\input psfig.sty
\input rotate

\magnification=1000
\tolerance=1000
\hsize=16truecm
\vsize=25truecm
\voffset=-1.5truecm
\newdimen\pageheight
\pageheight=\vsize
\tolerance=1000
\font\smain=cmss12 scaled\magstep0
\font\small=cmss10 scaled\magstep0
\font\sbig=cmss12 scaled\magstep2
\font\smedbf=cmssbx12 scaled\magstep1
\font\sbigbf=cmssbx12 scaled\magstep3

\def\oneplushalf{\baselineskip=18pt}

\parskip=12pt
\parindent=0pt
\oneplushalf

\def\leaderfill{\leaders\hbox to 1em{\hss.\hss}\hfill}
\def\ent #1 #2 #3{\baselineskip 12pt \line{#1 \hskip .6truepc #2 \leaderfill \quad #3}}
\def\presubent #1 #2{\baselineskip 12pt \line{\hskip 2.2truepc #1 \hskip .6truepc #2 \hfil}}
\def\subent #1 #2 #3{\baselineskip 12pt \line{\hskip 2.2truepc #1 \hskip .6truepc#2 \leaderfill \quad #3}}
\def\postsubent #1 #2{\baselineskip 12pt \line{\hskip 5.2truepc #1 \leaderfill \quad #2}}
\def\subsubent #1 #2 #3{\baselineskip 12pt \line{\hskip 5.2truepc #1 \hskip .6truepc #2 \leaderfill \quad #3}}

\nopagenumbers

{\smain
\hfill\vrule width2pt depth24.7cm
\vskip -\vsize
\vskip 3.0cm
\hfill {NCAR/TN--463+STR~~}
\vskip -14pt
\hfill {NCAR TECHNICAL NOTE~~}
\vskip 8pt
\hrule width7.0in height2pt
\vskip -4pt
\hfill {\small June 2004~~~}
\vskip 10pt

\sbigbf
{\noindent \sbigbf Scientific Description of the Sea Ice} 
\vskip 1pt
{\noindent \sbigbf Component in the Community Climate}
\vskip 1pt
{\noindent \sbigbf System Model, Version Three}
\vskip 12pt
\sbig
{\noindent B.~P.~Briegleb, C.~M.~Bitz, E.~C.~Hunke}
\vskip -10pt
{\noindent W.~H.~Lipscomb, M.~M.~Holland,}
\vskip -10pt
{\noindent J.~L.~Schramm, and R.~E.~Moritz}
\smain

\vskip 10cm

\vfill{
\hfill {CLIMATE AND GLOBAL DYNAMICS DIVISION~~~}
\vskip 8pt
\hrule width7.0in height2pt
\vskip -4pt
\hfill {NATIONAL CENTER FOR ATMOSPHERIC RESEARCH~~~}
\vskip -14pt
\hfill {BOULDER, COLORADO~~}
} }
\vfill\eject

\null

\vfill\eject

\parskip=6pt
\hsize=6.5truein
\tolerance=9000
\baselineskip=24pt
\footline={\centerline{\folio}}
\vsize=8.8truein
\voffset=0pt
\parindent=12pt

\font\main=cmr10 scaled\magstep1
\font\bf=cmbx10 scaled\magstep1
\font\it=cmsl10 scaled\magstep1

\main

\pageno=-1

\null

\vskip 2true cm
\centerline {\bf SCIENTIFIC DESCRIPTION OF THE SEA ICE}
\centerline {\bf COMPONENT IN THE COMMUNITY CLIMATE}
\centerline {\bf SYSTEM MODEL, VERSION THREE}

\vskip 2true cm
\centerline{ by Bruce P. Briegleb, Cecilia M. Bitz\footnote{$^1$}{University of Washington},
Elizabeth C. Hunke\footnote{$^2$}{T-3 Fluid Dynamics Group, Theoretical Division, Los 
Alamos NM 87545},}
\centerline{William H. Lipscomb$^2$, Marika M. Holland,}
\centerline{Julie L. Schramm, and Richard E. Moritz$^1$}
\vskip 0.75true cm
\centerline{National Center for Atmospheric Research}
\centerline{P.O. Box 3000, Boulder, CO 80307}

\vfill\eject

\null

\nopagenumbers

\vfill\eject

\footline={\ifnum\pageno>-8\line{\hfil\tenrm\folio\hfil}\fi}

\null

\pageno=-3

\vskip 2.5true cm

{\bf \centerline{TABLE OF CONTENTS}}
\vskip 12pt
\line{{Preface} \hskip .6truepc \leaderfill \quad {v}}
\vskip 10pt
\ent{1.} {Introduction} {1}
\vskip 8pt
\ent{2.} {Overview of the Community Sea Ice Model (CSIM)} {3}
\vskip 4pt
\subent{2.1} {State Variables} {3}
\vskip 4pt
\subent{2.2} {Fundamental Equations} {3}
\vskip 4pt
\subent{2.3} {Boundary Conditions and Solution} {6}
\vskip 4pt
\subent{2.4} {Summary} {9}
\vskip 8pt
\ent{3.} {Discretization} {10}
\vskip 4pt
\subent{3.1} {Time} {10}
\vskip 4pt
\subent{3.2} {Thickness} {10}
\vskip 4pt
\subent{3.3} {Vertical} {10}
\vskip 4pt
\subent{3.4} {Horizontal} {11}
\vskip 4pt
\subent{3.5} {Domain Decomposition} {12}
\vskip 8pt
\ent{4.} {Parameterizations and Numerical Approximations} {13}
\vskip 4pt
\subent{4.1} {Thickness Distribution} {13}
\vskip 4pt
\subent{4.2} {Thermal Properties} {14}
\vskip 4pt
\subent{4.3} {Input from the Coupler} {16}
\vskip 4pt
\subent{4.4} {Snow and Ice Albedo} {17}
\vskip 4pt
\subent{4.5} {Ice to Atmosphere Flux Exchange} {20}
\vskip 4pt
\subent{4.6} {Vertical Heat Conduction} {24}
\vskip 4pt
\subsubent{4.6.1} {Surface boundary conditions} {26}
\vskip 4pt
\subsubent{4.6.2} {Case I: Snow accumulation with no melting} {27}
\vskip 4pt
\subsubent{4.6.3} {Case II: Snow free with no melting} {29}
\vskip 4pt
\subsubent{4.6.4} {Case III: Snow accumulation with melting} {30}
\vskip 4pt
\subsubent{4.6.5} {Case IV: No snow with melting} {30}
\vskip 4pt
\subent{4.7} {Thickness Changes: Top, Bottom and Snow-to-Ice Conversion} {31}
\vskip 4pt
\subent{4.8} {Lateral Formation, Side and Bottom Melt Fluxes} {33}
\vskip 4pt
\subent{4.9} {Linear Remapping} {35}
\vskip 4pt
\subent{4.10} {Velocity} {38}

\vfill\eject
\null
\vskip 10pt
\vskip 2.5true cm

\subent{4.11} {Advection} {46}
\vskip 4pt
\subsubent{4.11.1} {Reconstructing area and tracer fields} {47}
\vskip 4pt
\subsubent{4.11.2} {Locating departure triangles} {48}
\vskip 4pt
\subsubent{4.11.3} {Integrating transports} {49}
\vskip 4pt
\subsubent{4.11.4} {Updating state variables} {50}
\vskip 4pt
\subent{4.12} {Mechanical Redistribution} {50}
\vskip 4pt
\subent{4.13} {Output to the Coupler} {56}
\vskip 8pt
\ent{5.} {Active Ice Only (AIO) Framework} {58}
\vskip 8pt
\ent{6.} {Output to History Files} {60}
\vskip 8pt
\ent{7.} {Summary} {64}
\vskip 8pt
\line{{Acknowledgments} \hskip .6truepc \leaderfill \quad {64}}
\vskip 8pt
\line{{List of Physical Constants} \hskip .6truepc \leaderfill \quad {65}}
\vskip 8pt
\line{{Appendix} \hskip .6truepc \leaderfill \quad {67}}
\vskip 8pt
\line{{References} \hskip .6truepc \leaderfill \quad {68}}
\vfill\eject

\null

\vfill

\parskip=8pt
\hsize=6.5truein
\tolerance=9000
\baselineskip=16pt
\footline={\centerline{\folio}}
\vsize=8.8truein
\voffset=30pt
\parindent=12pt

\pageno=-5

\null

\centerline {\bf Preface}
\noindent

Based on the Climate System Model Version 1 (CSM1) simulations 
(Boville and Gent, 1998; Weatherly et al., 1998), the CSM Polar 
Climate Working Group recommended a number of improvements 
to the sea ice component (the sea ice component is referred to as 
the Community Sea Ice Model, or CSIM). The most recent versions of the 
Community Climate System Model (CCSM), Versions Two and Three, include
improved versions of CSIM which satisfy all of those recommendations 
and include additional parameterizations, representing a major 
improvement over CSM1. The sea ice component of CCSM2, CSIM4, is
described in Briegleb et al. (2002) on the CCSM web page, and the 
enhanced sea ice component of CCSM3, CSIM5, is presented here.

CSIM5 consists of: elastic-viscous-plastic dynamics (Hunke and 
Dukowicz, 1997), which includes the effects of metric terms (Hunke and 
Dukowicz, 2002), energy conserving thermodynamics with a resolved 
vertical temperature profile and an explicit brine pocket parameterization 
(Bitz and Lipscomb, 1999), Lagrangian ice thickness distribution 
(Thorndike et al., 1975; Bitz et al., 2001), linear remapping for 
thickness space evolution (Lipscomb, 2001), mechanical redistribution 
due to rafting and ridging (Hibler 1980), ice strength computed from 
energetics (Rothrock, 1975), lateral and bottom melt processes (McPhee, 
1992), second order horizontal advection using remapping (Lipscomb and 
Hunke, 2004) and an albedo parameterization with implicit melt ponds. 
Five thickness categories adequately resolve the ice thickness distribution. 
Flux exchange with the atmosphere and ocean is evaluated over each 
thickness category and aggregated. 

CSIM5 uses two-dimensional domain decomposition and time split thermodynamics 
and dynamics for efficient parallel performance, and is vectorized for 
efficient vector performance. The code is written using standard parallel
and Fortran 90 constructs, and runs on several platforms. 

This document presents the details of the CSIM5 physical justifications, 
fundamental equations, parameterizations, numerical approximations/algorithms 
and run-time output.

\vfill\eject

\null

\nopagenumbers

\vfill\eject

\pageno=1

\footline={\centerline{\folio}}

\vskip 10pt
\noindent {\bf 1. Introduction}

     The Community Climate System Model (CCSM) is a coupled climate
model consisting of atmosphere, ocean, land, and sea ice components.
The model includes a coupler which passes fluxes and state variables
from one model component to another.  The model has been designed to 
support experiments that contribute to the understanding of past, 
present and future climates.  

     Because of the great annual range of insolation at the top of 
the atmosphere in high latitudes, the ocean surface there loses
heat to the atmosphere for months at a time, producing sub-freezing
surface temperature and a sea ice cover.  During the sunlit seasons,
sea ice persists because of its relatively high surface albedo and
its latent heat of formation.  During the winter, sea ice interposes
a layer of effective thermal insulation between the liquid ocean
and the atmosphere.  These effects keep the surface and lower atmosphere
over high latitude oceans much colder than they would be in the 
absence of sea ice, thereby affecting the horizontal gradient of 
air temperature and the general circulation of the atmosphere. 
When external forcing perturbs the earth's heat budget, the
sea ice cover may expand, contract, thicken or thin in such a 
way as to amplify the perturbation.  Therefore
it is important to include a realistic representation of sea ice
processes in models that would produce accurate simulations of 
past, present and future climate.

     The purpose of this document is to describe the sea ice component
of the current version of the CCSM.   This sea ice model is called
the Community Sea Ice Model (CSIM), and its current version
is CSIM5 which is being released as part of CCSM3. 
Emphasis is placed on relationships between the subcomponents
of CSIM5, and on aspects of the model of interest to users
who wish to perform simulations and diagnose the output. CSIM5
has evolved from earlier versions associated with the 
CCSM project.  CSIM1, which was released as part of the 
Climate System Model version 1 (CSM1) resolved one category of
sea ice thickness and included the following dependent 
variables:  ice concentration, ice thickness, snow depth, surface
temperature, ice temperature and ice velocity (Bettge, et al., 
1996).  The thermodynamic part of CSIM1 was based on Semtner (1976).
The dynamics were modeled as a cavitating fluid (Flato and 
Hibler, 1992).  Coupled experiments with CSM1 resulted in significant
biases in sea ice thickness and other climate variables
over the Arctic (Boville and Gent, 1998; Weatherly, et al., 1998).  

     The biases in CSIM1 resulted from the combined 
effects of the atmosphere, sea ice, ocean and land components
of the model, and were not uniquely attributed to errors in 
CSIM1.  Nevertheless, the biases called attention to simplifications
and approximations in CSIM1 at a time when considerably 
more elaborate representations of ice dynamics and thermodynamics 
had been developed for regional ice-ocean models.  
Shortly after the release of CSIM1 the CCSM Polar Climate 
Working Group (PCWG) identified aspects of the sea ice model
that should be given high priority in development efforts:
(1) A plastic rheology with an elliptical yield curve;
(2) Enhanced sea ice thermodynamics; (3) A multi-category
ice thickness distribution; (4) Elimination of spurious 
ice convergence at the North Pole; (5) Compatibility of 
the grids used to model the sea ice and the ocean; 
(6) Efficient parallelization of the model code; 
(7) Development of an active ice-only version framework
for testing the model without full coupling to the 
atmosphere and ocean model components.  These seven
high-priority objectives were achieved with the 
release in 2002 of CSIM4 as a component of CCSM2.
A description of CSIM4 is available on the CCSM website
(Briegleb, et al., 2002).

      After the release of CSIM4 the PCWG decided to give 
high priority to developing an efficient vector version of the
code, and to enhance some of the numerical algorithms in the 
model, resulting in CSIM5.  In CSIM5 the second order horizontal
advection scheme (Smolarkiewicz, 1984) has been replaced by
a more accurate and numerically efficient remapping 
scheme (Lipscomb and Hunke, 2004).  The dynamic boundary 
condition for marginal ice zones has been improved (Connolley 
et al. 2004).  Salinity is now exchanged explicitly between 
the ocean and the sea ice model components.  Finally, the 
albedo of snow and ice surfaces have been adjusted to 
accomodate changes in the simulation of the polar atmosphere
by CCSM3.  

     This document is organized as follows.  Section 2 provides
an overview of CSIM5 by listing the state variables, fundamental
equations, boundary conditions and a summary of the model physics.
Section 3 presents the time and space discretizations of the 
fundamental equations.  Section 4 is the longest section,
providing detailed descriptions of the parameterizations and
numerical approximations.  Sections 5, 6, and 7 present the 
active ice-only framework, the output to history files,
and a summary of CSIM5 respectively.  In addition to the supported,
default physics of the CSIM5 there are alternative physics
options available, including developmental versions of the
model.   These options are summarized in the Appendix. 

    Additional information on how to obtain and run CSIM5 is available 
in the "CSIM User's Guide". For details on source code structure and 
its modification, see the "CSIM Code Reference". Both documents can be 
found on the CCSM web page under models (http://www.ccsm.ucar.edu/models),
as well as the present document.
\vfill\eject

\vskip 10pt
\noindent {\bf 2. Overview of the Community Sea Ice Model (CSIM)}

\vskip 8pt
\noindent {\bf 2.1 State Variables}

The {\bf{state variables}} for the sea ice model are listed in Table 1. Where possible, 
we use conserved quantities as state variables. 

\vskip 15pt
\vbox{
\hang \noindent {{\it Table 1.} 
State Variables for the sea ice model. Subscript $n$, $\{n=0,1,2,..N\}$ refers 
to the $n^{th}$ ice thickness category ($n=0$ is open water), 
where $N$ is the total number of categories. 
For CSIM5, $N=5$. Subscript $l$, $\{l=1,..L\}$ refers to vertical level, with $L=4$. 
For each category, ice thickness 
($h_n=\Sigma_{l=1}^L V_{nl} / A_n$) lies within constant category thickness limits.
Ice velocity $\bf u$ and the associated stress tensor $\sigma$ (components 
$i$=1,2; $j$=1,2) are not resolved across the ice thickness distribution. }
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Symbol & Description \cr
\noalign{\hrule height 1pt}
         $A_n$                  & Sea ice area (fraction from 0 to 1)          \cr
         $V_{nl}$               & Sea ice volume  (m$^3$ m$^{-2}$)             \cr
         $E_{nl}$               & Sea ice internal energy (J m$^{-2}$)         \cr
         $V_{sn}$               & Snow volume (m$^3$ m$^{-2}$)                 \cr
         $T_{sn}$               & Surface temperature of snow/ice ($^\circ$C)  \cr
         $$\bf u$$              & Sea ice velocity (m s$^{-1}$)                \cr
         $\sigma_{ij}$          & Stress tensor components (N m$^{-1}$)        \cr
\noalign{\hrule height 1pt}
}}}
\vskip 10pt

\vskip 8pt
\noindent {\bf 2.2 Fundamental Equations}

The {\bf{fundamental equations}} determine the spatial and temporal evolution
of the state variables. We first give a rationale for using an ice
thickness distribution before describing these equations.

Many properties of sea ice depend on ice thickness (Thorndike et al.,1975);
for example, ice compressive strength, growth rate, surface temperature, 
turbulent and radiative flux exchange with the atmosphere. Two contrasting 
phenomena alter the distribution of ice thickness on a yearlong average: 
accretion and ablation against lead opening and ridging of ice. This 
competition was elegantly described by Thorndike et al. (1975): ``thermodynamics 
seeks the mean and dynamics the extremes''. The evolution of the thickness 
distribution is the historical integral of these two continuous processes.

Formally, the thickness distribution is described by the distribution
function $g(h,{\bf{x}},t)$, where $h$ is ice thickness (henceforth we 
suppress the explicit space and time dependence). $g(h)dh$ is the 
fraction of area covered by ice of thickness $h$ to $h+dh$, normalized 
by $\int_0^{\infty} g(h) dh = 1$, the conservation of total area. 
The aggregate ice fraction is $A = \int_{0^+}^{\infty} g(h) dh$,
while the open water fraction is $A_0 = g(h=0) = 1 - A$. The cumulative 
distribution function is $G(h) = \int_0^h g(h) dh$. The average of 
a quantity $F$ that depends on ice thickness is referred to as the 
aggregate ${\bar F} = {1 \over A}\int_0^{\infty} F(h)g(h) dh$.

The evolution of $g$ is governed by the distribution equation
$$
         {{\partial g} \over {\partial{t}}} = 
- {{\partial} \over {\partial h}} ({\dot{h}}g)
+ L(h,g) 
- \bigtriangledown\cdot ({\bf u} g)  
+ R(h,g,{\bf u})   \eqno(1)
$$
where $\dot h$ is the rate of change in ice thickness due to vertical thermodynamic 
processes, $-{\partial \over \partial h} ({\dot{h}}g)$ is the change in distribution
due to thickness space transport, 
$L(h,g)$ is the change in distribution due to lateral 
melt/formation processes, $- \bigtriangledown\cdot ({\bf u} g)$ is the change
in distribution due to horizontal advection ($\bigtriangledown$ is the horizontal 
gradient operator and $\bf u$ is the velocity field over the 
thickness distribution) and $R(h,g,{\bf u})$ ($\psi$ in Thorndike et al.,1975)
is a redistribution function due to rafting and ridging processes.

To solve Eq. 1, a discrete set of N ice categories is assumed, delimited 
by the thicknesses $h_n^\ast,\{n=0,1,2...N\}$ for which $h_0^\ast=0$. Thus,
Eq. 1 is integrated over the thickness limits for each category, 
resulting in a discrete set of N equations to be solved for the ice fraction
in each category n:
$$
A_n = \int_{h_{n-1}^\ast}^{h_n^\ast} g(h) dh,             \eqno(2)
$$
where the total (aggregate) ice fraction $A = \sum_{n=1}^N A_n$. The first moment 
of the distribution function is the ice volume $V_n = \int_{h_{n-1}^\ast}^{h_n^\ast} 
hg(h) dh$, and any function $F$ which is linear in the ice thickness (i.e. $F = F_0 + 
F_1h$) results in $F_n = F_0 A_n + F_1 V_n$.

One way to solve Eq. 1 is to assume thickness is distributed uniformly within each 
category (Hibler, 1980; Flato and Hibler, 1995), resulting in Eulerian advection in 
thickness space due to sea ice growth and melt processes. Such advection is very 
diffusive unless a large number of categories are employed. A second way to solve 
Eq. 1 is to assume ice can vary in thickness within each category (Thorndike 
et al., 1975). This Lagrangian method is free from the diffusion of the Eulerian 
thickness advection, allowing for a smaller number of categories, as well as the 
resolution of the vertical temperature profile in snow and ice (Bitz et al., 2001). 
It is used for the present sea ice model, except for the linear remapping which 
uses a combination of both Eulerian and Langrangian methods.

Thus, the fundamental equations for the present sea ice model start with the
discrete form of Eq. 1 for the ice fractions $A_n$, the first moment 
equations for the ice volume $V_n$, corresponding equations for the snow volume
$V_{sn}$, equations for the vertically varying ice internal energy (from which
the vertical temperature profile and heat transfer are evaluated), equations
for the surface temperature required for the vertical heat transfer solution,
and finally dynamic equations for the ice velocity $\bf{u}$
needed to evaluate horizontal advection and the ridging terms in the 
distribution equations.

The fundamental equations are as follows. For the category sea ice fraction
and volume:
$$ 
         {\partial A_n \over \partial t} = S_{TAn} - \bigtriangledown\cdot 
         (\bf{u}A_n) + S_{MAn}~~~(n=1,2,...5)                     \eqno(3)
$$
$$
         {\partial V_n \over \partial t} = S_{TVn} - \bigtriangledown\cdot 
         (\bf{u}V_n) + S_{MVn}~~~(n=1,2,...5)                     \eqno(4)
$$
where terms $S_T$ denote sources/sinks due to thermodynamic processes and 
thickness space transport, while terms $S_M$ denote sources/sinks due to 
mechanical redistribution. The sea ice thickness $h_n$ is derived from 
the fraction and volume as $h_n=V_n/A_n$.

To resolve vertical atmosphere/ocean heat exchange, and account as well for internal 
heat within the ice, the ice internal energy $E_n$ (vertically varying) is governed 
by the conservation equation:
$$
         {\partial E_n \over \partial t} = S_{TEn} - \bigtriangledown\cdot 
         (\bf{u}E_n) + S_{MEn}~~~(n=1,2,...5)                      \eqno(5)
$$
The internal sea ice energy $E_n$ is proportional to the ice volume, $E_n = q_nV_n$, 
where the proportionality function $q_n$ (termed the energy of melting, or enthalpy) 
is the internal
energy per unit volume. The effects of brine pockets are represented explicitly through 
the temperature $T_n$ and salinity $S_n$ dependent energy of melting $q_n=q_n(T_n,S_n)$. 
The vertical temperature profile is inferred by solving for $T_n$ in $q_n(T_n,S_n)=E_n/V_n$ 
over an ice thickness $h_n=V_n/A_n$, using a prescribed salinity profile $S_{n}$.

The conservation equation for snow volume $V_{sn}$ is:
$$
         {\partial V_{sn} \over \partial t} = S_{TVsn} - \bigtriangledown\cdot 
         (\bf{u}V_{sn}) + S_{MVsn}~~~(n=1,2,...5)                   \eqno(6)
$$
Snow thickness is derived from $h_{sn}=V_{sn}/A_n$. Snow energy per unit volume 
is $E_{sn} = q_s V_{sn}$, where the energy of melting of snow $q_s$ is constant.

For each category, the heat equation governing vertical heat transfer over time 
interval $t$ to $t'$ corresponding to temperatures $T_n$ and $T'_n$ respectively, 
allowing for temperature and salinity dependent heat capacity $c_i$, thermal 
conduction and internal absorption of penetrating solar radiation, is given by:
$$
  \int^{T'_n}_{T_n} \rho_i c_i dT_n = 
  \int^{t'}_t \left( { \partial  \over \partial z} k { \partial T_n 
  \over \partial z} + Q_{SW} \right) dt~~~(n=1,2,...5)  \eqno(7)
$$
where sea ice is assumed to have a constant density $\rho_i$, $z$ is the vertical 
coordinate within the sea ice, $Q_{SW}$ is the absorbed shortwave flux, and the 
thermal conductivity $k$ is that for either snow or ice. Modifications to the 
temperature profile resulting from heat transfer change the ice internal energy 
according to $E_n=q_n(T_n,S_n)V_n$.

The surface boundary conditions for the vertical heat transfer solution require 
surface temperature $T_{sn}$ to satisfy the conservation equation
$$
         {\partial A T_{sn} \over \partial t} = S_{TTsn} - \bigtriangledown\cdot 
         (\bf{u} A_n T_{sn}) + S_{MTsn}~~~(n=1,2,...5)            \eqno(8)
$$
Evaluation of the thermodynamic source term $S_{TTsn}$ at the surface requires
calculation of surface snow and ice albedo, which are diagnostic functions of
solar spectral interval, snow and ice thickness, and surface temperature.

For momentum conservation, sea ice is assumed to be a two-dimenstional
continuum. The sea ice velocity $\bf{u}$ and stress tensor $\sigma_{ij}$ are 
considered (along with related dynamic quantities) to be representative of the 
entire ice thickness distribution. Their governing equations are:
$$
        {\bar m} {\partial \bf{u} \over \partial t} = 
          -{\bar m}f\bf{k}\times{\bf{u}} + \tau_a + \tau_o + 
           {\bar m}g_e\nabla{H_o} + \nabla\cdot\sigma   \eqno(9) 
$$
where ${\bar m} = \rho_s {\sum_{n=1}^N}V_{sn} + \rho_i {\sum_{n=1}^N}V_n $,
the non-linear $\bf{u}$ advection terms are ignored as they are 
negligibly small when the equations are scaled, $f$ is the Coriolis parameter, 
$\bf{k}$ is the local vertical unit vector,
$\tau_a$ and $\tau_o$ are air and water stresses respectively, $g_e$ is 
the gravitational acceleration, $H_o$ is the sea surface height 
and $\nabla\cdot\sigma$ is the force per unit area due to internal 
ice stress, where $\sigma$ is the stress tensor. 
The stress tensor equations are:
$$
      {\partial \sigma_{ij} \over \partial t} + 
      {e^2 \over 2T_{ew}} \sigma_{ij} + 
      {1 - e^2 \over 4T_{ew}} \sigma_{kk} \delta_{ij} = 
        {P \over 2T_{ew}\Delta^\prime} {\dot\epsilon_{ij}} 
       - {P \over 4T_{ew}} \delta_{ij}~~~(i,j=1,2)               \eqno(10) 
$$
where $(i,j = 1,2)$ refer to the four components of the stress tensor, $e$ is a
constant ratio of major to minor axes of the elliptical yield curve, $T_{ew}$ is a 
damping time scale for elastic waves, $\delta_{ij}$ is the Kronecker delta, $P$ 
is the ice compressive strength (or mechanical pressure, a function of the thickness 
distribution), ${\dot\epsilon_{ij}}$ is the rate of strain tensor, in turn a function 
of velocity gradients, and $\Delta^\prime$ is a function of the rate of strain tensor.

\vskip 8pt
\noindent {\bf 2.3 Boundary Conditions and Solution}

{\bf{Boundary conditions}} are represented vertically by atmospheric/oceanic 
forcing from the coupler, consisting of states and interfacial fluxes summarized 
in Table 2, and horizontally by no-slip $\bf{u} \rightarrow 0$ along coastlines 
and $\bf{u} \rightarrow \bf{u_o}$ (ocean surface current) on the open ocean edge. 
The fundamental equations are solved subject to these boundary conditions, and
selected states along with atmosphere and ocean fluxes are returned to the coupler
as listed in Table 3.

The atmospheric states are used, along with ice surface temperature and roughness, 
to compute the bulk sensible/latent heat fluxes $F_{SH}, F_{LH}$ and surface stress
components $\tau_{ax},\tau_{ay}$. 
The downwelling shortwave flux is partitioned directionally and spectrally into four 
components such that the total downwelling shortwave flux is 
$F_{SWDN} = F_{SWvsdr} + F_{SWvsdf} + F_{SWnidr} + F_{SWnidf}$, while the total 
downwelling longwave flux is $F_{LWDN}$.

\vskip 15pt
\vbox{
\hang \noindent {{\it Table 2.} 
State Variables and Fluxes Received by Sea Ice Model from the Coupler}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
\noalign{\hrule height 1pt}
 Symbol & Description & Units \cr
\noalign{\hrule height 1pt}
                        & Atmospheric States                              &                \cr
\noalign{\hrule height 1pt}
         $z_a$          &  reference height                               &  m             \cr
         $u_a$          &  x direction wind speed at $z_a$                &  m s$^{-1}$    \cr
         $v_a$          &  y direction wind speed at $z_a$                &  m s$^{-1}$    \cr
         $\theta_a$     &  potential temperature at $z_a$                 &  K             \cr
         $T_a$          &  air temperature at $z_a$                       &  K             \cr
         $q_a$          &  specific humidity at $z_a$                     &  kg kg$^{-1}$  \cr
         $\rho_a$       &  air density at $z_a$                           &  kg m$^{-3}$   \cr 
\noalign{\hrule height 1pt}
                        & Atmosphere $\Rightarrow$ Sea Ice Fluxes (+ downwards) &  \cr
\noalign{\hrule height 1pt}
         $F_{SWvsdr}$   &  direct, visible downwelling shortwave          &  W m$^{-2}$  \cr
         $F_{SWvsdf}$   &  diffuse, visible downwelling shortwave         &  W m$^{-2}$  \cr
         $F_{SWnidr}$   &  direct, near infrared downwelling shortwave    &  W m$^{-2}$  \cr
         $F_{SWnidf}$   &  diffuse, near infrared downwelling shortwave   &  W m$^{-2}$  \cr
         $F_{LWDN}$     &  downwelling longwave                           &  W m$^{-2}$  \cr
         $F_{RN}$       &  water flux due to rain                         &  kg m$^{-2}$ s$^{-1}$  \cr
         $F_{SNW}$      &  water flux due to snow (liquid equivalent)     &  kg m$^{-2}$ s$^{-1}$  \cr 
\noalign{\hrule height 1pt}
                        & Ocean States                                    &  \cr  
\noalign{\hrule height 1pt}
         $T_{ocn}$      &  sea surface temperature                        &  K            \cr
         $S_{ocn}$      &  sea surface salinity                           &  ppt          \cr
         $u_o$          &  x direction ocean surface current              &  m s$^{-1}$         \cr
         $v_o$          &  y direction ocean surface current              &  m s$^{-1}$         \cr
   $(\nabla H_o)_x$     &  x direction sea surface slope                  &  m m$^{-1}$      \cr
   $(\nabla H_o)_y$     &  y direction sea surface slope                  &  m m$^{-1}$      \cr 
\noalign{\hrule height 1pt}
                        & Ocean $\Rightarrow$ Sea Ice Fluxes (+ downwards) &  \cr  
\noalign{\hrule height 1pt}
         $F_{Qoi}$      &  freezing/melting potential                     &  W m$^{-2}$       \cr 
\noalign{\hrule height 1pt}
}}}

Based on the ice state, snow/ice directional and spectral albedos $\alpha_{vsdr},
\alpha_{vsdf},\alpha_{nidr},\alpha_{nidf}$ are evaluated, and used to compute the
total absorbed shortwave in the ice $F_{SW}$. Of this total a portion $I_{SW}$ 
penetrates below the surface and is either internally absorbed in the ice 
$Q_{SW}$ or penetrates to the ocean below as $F_{SWo}$. The net longwave flux at 
the surface $F_{LW}$ is the difference between the upwelling longwave $F_{LWUP}$
and the absorbed downwelling longwave. The upwelling longwave is given by the
the surface emission $\varepsilon \sigma_{sb} T_s^4$ (where $\varepsilon$ is the 
snow/ice longwave emissivity and $\sigma_{sb}$ is the Stefan-Boltzmann constant) 
and the reflection of downwelling longwave $(1-\varepsilon) F_{LWDN}$, while the 
absorbed downwelling longwave is $\varepsilon F_{LWDN}$.

Rain $F_{RN}$ is assumed to run off directly into the ocean, and thus 
contributes to the ocean-ice water flux $F_{Wo}$. Snow $F_{SNW}$ 
is used to compute snow accumulation $dh_s/dt = F_{SNW} / \rho_s$. 
The latent heat flux $F_{LH}$ is associated with an evaporative
water flux to the atmosphere $F_{EVAP}$, and an associated salt
flux $S_i \rho dh/dt$ with the ocean, where $h$ is ice thickness and $S_i$ 
is a constant sea ice reference salinity, which contributes to the
ocean-ice salt flux $F_{So}$. This salt exchange with the
ocean for sublimation/condensation is required to maintain a constant 
sea ice reference salinity.

The top boundary condition at surface temperature $T_s$ is
$F_{TOP}(T_s) = F_{SW}-I_{SW} + F_{LW} + F_{SH} + F_{LH} + k dT/dz$, while 
the lowest ice interface in contact with the ocean is at ocean freezing 
temperature. With these boundary conditions and internal shortwave heating,
the heat equation can be solved. If $F_{TOP}(T_{melt}) > 0$, where $T_{melt}$ 
is the snow/ice melting temperature, then snow/ice melt is computed by 
$F_{TOP}(T_s=T_{melt}) = q dh/dt$ as appropriate for either snow (if present) 
or ice ($h$ is the thickness of snow or ice). Snow and ice melt is assumed to run 
off directly into the ocean, contributing to the ocean-ice water flux $F_{Wo}$ 
by the amount $\rho dh/dt$ where $\rho h$, are the density,thickness of snow or 
ice, respectively. If sea ice melts, then an additional salt flux of 
$S_i \rho dh/dt$ is exchanged with the ocean.

Ice formation occurs by three processes. Although these processes are distinguished
in formation, no distinction is made between ice types. If the freezing/melting
potential ($F_{Qio}$ in Table 2.) is such that heat is required by the ocean to maintain 
the freezing temperature ($F_{Qoi}>0$), then {\bf{frazil ice}} formation occurs, at a 
rate $dV_f/dt = F_{Qoi} / \rho_i q_f$, where $q_f$ is a heat of melting assuming the
ice forms at $T=0^\circ C$ and $S=0$. Frazil ice formation has an implied salt 
flux to the ice, sufficient to establish $S_i$ as the mean salinity of the newly 
formed sea ice. If the freezing/melting potential indicates ocean heat is available 
to melt ice $F_{Qoi}<0$, this heat is partitioned between lateral $F_{SID}$ and bottom 
$F_{BOT}$ heat fluxes according to the fraction of absorbed solar energy near the 
surface and in deeper water. The sea surface temperature $T_{ocn}$ is used to compute heat 
fluxes $F_{SID}$ and $F_{BOT}$ (sea surface salinity $S_{ocn}$ is not used). The heat flux 
for lateral and bottom melting is $F_{Qio}$, with associated water and salt fluxes 
$\rho_i dh/dt$ and $S_i\rho dh/dt$ to the ocean, respectively.

The bottom boundary condition is
$F_{BOT} - k dT/dz = q dh/dt$. If bottom ice formation occurs (i.e. $dh/dt > 0$),
this ice is termed {\bf{congelation ice}}. If sufficient snow $h_s$ overlies 
ice, the snow-ice interface can be depressed below sea level. Snow below sea level 
is converted into ice conserving mass and energy, and is termed {\bf{snow-ice}}. For 
both of these processes, a salt flux of $S_i \rho dh/dt$ is exchanged with the ocean.

The ocean surface currents $(u_o,v_o)$ and ice velocity $(u,v)$ are used
to compute ocean/ice stress $(\tau_{ox},\tau_{oy})$. The tilt stress is computed 
from the gradient of the sea surface height ($(\nabla H_o)_x$,$(\nabla H_o)_y$).
With these stresses, sea ice velocity and ridging are computed.

\vskip 15pt
\vbox{
\hang \noindent {{\it Table 3.} 
State Variables and Fluxes Sent from Sea Ice Model to the Coupler}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
\noalign{\hrule height 1pt}
 Symbol & Description & Units \cr
\noalign{\hrule height 1pt}
                        & Sea Ice States                           &       \cr
\noalign{\hrule height 1pt}
         $A$            &  ice area                                & fraction (0 to 1) \cr
         $T_s$          &  surface temperature                     &  K    \cr 
         $\alpha_{vsdr}$&  albedo (visible, direct)                & fraction (0 to 1) \cr
         $\alpha_{vsdf}$&  albedo (visible, diffuse)               & fraction (0 to 1) \cr
         $\alpha_{nidr}$ &  albedo (near infrared, direct)         & fraction (0 to 1) \cr
         $\alpha_{nidf}$ &  albedo (near infrared, diffuse)        & fraction (0 to 1) \cr
\noalign{\hrule height 1pt}
                        & Sea Ice $\Rightarrow$ Atmosphere Fluxes (+ downwards) &  \cr
\noalign{\hrule height 1pt}
         $F_{LH}$       &  latent heat flux                        &  W m$^{-2}$    \cr 
         $F_{SH}$       &  sensible heat flux                      &  W m$^{-2}$    \cr 
         $F_{LWUP}$     &  upwelling longwave                      &  W m$^{-2}$    \cr 
         $F_{EVAP}$     &  evaporated water                        &  kg m$^{-2}$ s$^{-1}$ \cr 
         $\tau_{ax}$    &  x direction atmosphere-ice stress       &  N m$^{-2}$    \cr 
         $\tau_{ay}$    &  y direction atmosphere-ice stress       &  N m$^{-2}$    \cr  
\noalign{\hrule height 1pt}
                        & Sea Ice $\Rightarrow$ Ocean Fluxes (+ downwards) &  \cr  
\noalign{\hrule height 1pt}
         $F_{SWo}$      &  shortwave transmitted to ocean          &  W m$^{-2}$    \cr 
         $F_{Qio}$      &  heat flux to ocean                      &  W m$^{-2}$    \cr 
         $F_{Wo}$       &  water flux                              &  kg m$^{-2}$ s$^{-1}$ \cr 
         $F_{So}$       &  salt flux                               &  kg m$^{-2}$ s$^{-1}$ \cr 
         $\tau_{ox}$    &  x direction ice-ocean stress            &  N m$^{-2}$    \cr 
         $\tau_{oy}$    &  y direction ice-ocean stress            &  N m$^{-2}$    \cr  
\noalign{\hrule height 1pt}
                        & Diagnostic Fields                        &  \cr  
\noalign{\hrule height 1pt}
         $T_{ref}$      &  reference temperature (2 m)             &  K             \cr 
         $Q_{ref}$      &  reference specific humidity (2 m)       &  kg/kg         \cr 
         $F_{SW}$       &  ice/ocean absorbed shortwave flux       &  W m$^{-2}$    \cr  
\noalign{\hrule height 1pt}
}}}

\vskip 8pt
\noindent {\bf 2.4 Summary}

Subject to initial and boundary conditions, Eqs. 3-10 constitute the 
fundamental equations for the sea ice model. The next two sections present 
the discretizations, parameterizations and numerical approximations used 
to solve the fundamental equations. 

\vfill
\eject
\noindent {\bf 3. Discretization}

\vskip 8pt
\noindent {\bf 3.1 Time}

The time stepping loop in the sea ice model is split into two intervals for 
improved CCSM computational performance. The physical time step, $\Delta t$, 
represents the time interval from the beginning to the end of the time step, 
whereas the coupling time step is displaced from this and denotes the time 
interval at which information is sent to the coupler. The order of calculations 
in the sea ice model differs somewhat from the order presented in Section 4. 
For more information on the time stepping, see the CSIM Code Reference.

During the first half of the physical time step, forcing fields in Table 2 are 
received from the coupler and the vertical thermodynamics are calculated. At 
this point, the ice surface state and ice-atmosphere fluxes in Table 3 are 
returned to the coupler. When the atmospheric model receives the ice-atmosphere 
fluxes, it can run in parallel with the ice model during the second half of the 
physical time step.

During the second half of the physical time step, lateral thermodynamics, thickness 
space transport, dynamics and physical space transport, mechanical and thermodynamic 
redistribution and albedo calculations are done. Note that since the lateral 
thermodynamics (including exchange with underlying ocean) and the albedo calculation 
follow the send to the coupler, the ice-ocean fluxes and snow/ice albedos 
are offset by one time step with respect to the ice states and ice-atmosphere fluxes. 
The albedos are computed last to ensure that the atmospheric computation  of the 
atmosphere-ice radiative fluxes received on the next time step use the {\bf{same}} 
albedos as are used for the ice vertical thermodynamic calculation, thus ensuring
energy conservation. 

\vskip 8pt
\noindent {\bf 3.2 Thickness}

The thickness distribution function $g(h)$ is integrated over $N$ discrete thickness 
ranges or categories. Presently there are $N=5$ ice thickness categories in the
standard model, but arbitrary $N$ is allowed. These categories are described in the
next section. The state variables (listed in Table 1) of
sea ice concentration, volume, energy, snow volume and surface temperature are 
discretized into $n=1,2,...N$ categories. Henceforth, a subscript $n$ will refer 
to the $n^{th}$ thickness category.

\vskip 8pt
\noindent {\bf 3.3 Vertical}

To compute vertical heat conduction through ice, ice thickness is divided into four
vertical layers. This requires sea ice internal energy $E$ to vary in the 
vertical over four evenly spaced layers in each thickness category. Temperatures are computed 
from $E$ using the energy of melting and the ice volume in each layer. Internal temperatures 
are centered within each layer, while conductivities and energy fluxes are represented at 
layer interfaces. Temperature boundary conditions at the surface and base of ice are taken 
at the top and bottom interfaces respectively. 

\vskip 8pt
\noindent {\bf 3.4 Horizontal}

The horizontal grids used for the sea ice model are displaced pole grids, 
in which the South Pole is typically located at the geographic South Pole, but the North
Pole may be located in any northern hemisphere land mass. The CCSM uses a grid with
the North Pole in central Greenland. The grids are orthogonal curvilinear, so that vectors 
parallel to increasing longitude and latitude coordinates are perpendicular to one 
another. Two available resolutions are the standard gx1 ($320$ longitudes x $384$ latitudes,
$\sim 1.1^\circ \times 0.94^\circ$), and a coarse gx3 ($100$ longitudes x $116$ latitudes,
$\sim 3.6^\circ \times 1.6^\circ$). The grids south of the equator are regular 
latitude/longitude spherical coordinates.

Spatial discretization is that of a B-grid. All state 
variables except ice velocity and stress tensor components are taken at grid box 
mid-points (a tracer grid termed the T-grid), while velocities are defined at grid box corners 
(a velocity grid termed the U-grid). The stress tensor, rates of strain and viscosities 
are defined bilinearly across each grid cell using the values at the corners. This
discretization tends to avoid decoupling problems associated with the B-grid. 

Grid information is taken from a grid data file read in by the ice model at
initialization. The fields on the U-grid are latitude, longitude and angle 
the grid makes with a geographic latitude line (ULAT, ULON and ANGLE respectively), 
and fields on the T-grid are the land/ocean mask and T cell widths on north and
east sides (tmask, HTN and HTE respectively). Land points on the T-grid are designated
by tmask, which is either 0 (land) or 1 (ocean). On the U-grid umask designates 
ocean points where ice velocities are possible, and is 0 if any of the 
four surrounding tmask points are 0 (i.e. land), and 1 otherwise. In other
words, umask is zero for all coastal and land points. The T-grid longitudinal and 
latitudinal widths through cell centers are dxt and dyt respectively, while 
dxu and dyu are the analogous widths for the U-grid. ANGLET ($\chi^t$) is ANGLE 
($\chi^u$) interpolated to the T-grid. Specifically:
$$
\eqalign{
      dxt_{ij}    &= {1 \over 2} (HTN_{ij}+HTN_{ij-1})    \cr
      dyt_{ij}    &= {1 \over 2} (HTE_{ij}+HTE_{i-1j})    \cr
      dxu_{ij}    &= {1 \over 2} (HTN_{ij}+HTN_{i+1j})    \cr
      dyu_{ij}    &= {1 \over 2} (HTE_{ij}+HTE_{ij+1})    \cr
      A^t_{ij}    &= dxt_{ij}~dyt_{ij}                    \cr
      A^u_{ij}    &= {1 \over 4} (A^t_{ij}+A^t_{i+1j}+A^t_{ij+1}+A^t_{i+1j+1}) \cr
      \chi^t_{ij} &= {1 \over 4} (\chi^u_{ij}+\chi^u_{i-1j}+\chi^u_{ij-1}+\chi^u_{i-1j-1}) \cr
} \eqno(13)
$$
where the $ij$ are the grid longitude and latitude indices respectively,
$A^t_{ij}$,$A^u_{ij}$ are the grid box areas on the T-grid and U-grid respectively, 
and for the $\chi^t$ calculation, $\chi^u$ values are adjusted if 
any differ by more than $180^\circ$.

\vskip 8pt
\noindent {\bf 3.5 Domain Decomposition}

The horizontal computational grid is domain decomposed in two dimensions for 
parallelization. The global domain of dimensions $imt_{global} \times jmt_{global}$ 
(for example, the gx1 grid has $imt_{global} = 320$ longitude points and $jmt_{global} 
= 384$ latitude points) is divided into integral NX longitude by NY latitude subdomains 
of dimensions ($imt_{local} = imt_{global}/NX + 2n_{ghost} + 1$) $\times$ ($jmt_{local} = 
jmt_{global}/NY + 2n_{ghost} + 1 $), where $imt_{global}/NX$, $jmt_{global}/NY$ must 
be integers. Each subdomain has a physical portion indexed as $[ilo:ihi,jlo:jhi]$ with 
$n_{ghost}$ boundary cells outside. Periodic boundary conditions are applied, with 
boundary routines performing communications between subdomains when running parallel. 
Global scatter and gather routines distribute information from the global domain to 
the subdomains and back.

We note that since the thermodynamic calculations involve one grid point at
a time, a purely thermodynamic model integration is independent of the domain
decompostion (i.e. the exact values of NX,NY), while the dynamic calculation 
depends upon domain boundary conditions. Hence an integration with active dynamics
is dependent upon the exact values of NX and NY.

\vfill
\eject
\noindent {\bf 4. Parameterizations and Numerical Approximations}

Sections 2 and 3 introduced the sea ice model state variables, fundamental 
equations, boundary conditions, solutions and discretizations. In the present 
section we elaborate both on the parameterizations necessary to represent 
various forcing terms in the fundamental equations and the details of the numerical 
solutions. For many of the subsections to follow, the processes are described 
(sometimes implicitly) for a particular thickness category $n$. For exchange with 
the atmosphere and the ocean however, only aggregate quantities (i.e. those summed 
over the thickness distribution) are used. 

\vskip 8pt
\noindent {\bf 4.1 Thickness Distribution}

The number of ice thickness categories $N$ (see Table 1) used in the ice model 
results from a trade off between the desire to resolve thin ice that is important 
in ocean-atmosphere heat exchange and feedback processes against computational
cost. Bitz et al. (2001) showed that five thickness categories with adequate thin
ice resolution are sufficient to represent the first order effects of an ITD.
We therefore chose $N=5$, with thickness boundaries given in Table 4,
based on the category limit formula of Lipscomb (2001).
While $N=5$ is a convenient value for climate modeling, it is not hardwired 
in the code; other values can be used, including $N=1$. Also, while the 
category $1$ lower limit is $0$, for thermodyanmic stability the minimum ice
thickness in category 1 is $h_{min}$. If ice thickness in this category falls
below $h_{min}$, it is reshaped using $h_{min} A_1^\prime = h_1 A_1$, where
$A_1^\prime$ is the adjusted category 1 ice fraction.

\vskip 15pt
\vbox{
\hang \noindent {{\it Table 4.} 
Ice Thickness Distribution (N=5)}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 $n$ & Range ($m$) \cr
\noalign{\hrule height 1pt}
         0      &  0            \cr
         1      & $0^+$ - 0.65  \cr
         2      & 0.65 - 1.39   \cr
         3      & 1.39 - 2.47   \cr
         4      & 2.47 - 4.60   \cr
         5      & $>$ 4.60      \cr
\noalign{\hrule height 1pt}
}}}
\vskip 10pt

In the Lagrangian method of ice thickness distribution used here, 
there is no process (thermodynamic or dynamic) that absolutely 
prevents ice from outgrowing its thickness range for a given 
category. While it is true that the incremental linear remapping
used to evaluate thickness space transport in most cases prevents
ice from outgrowing its category thickness limits, an occasional 
adjustment is still necessary, which is termed ``thermodynamic 
redistribution''. This process contributes to the thermodynamic 
source terms $S_T$ in the conservation equations (see Section 1).

Any category of ice which outgrows its upper thickness limit
is combined with the next thickness category. The combination
is done preserving ice area, volume, energy and snow volume. 
Any category of ice which melts below its lower thickness
limit, will be combined with the next lower category in a
manner similar to outgrowth just described, except for the
thinnest ice category. For ice in category 1 that melts below
its lower limit, the ice is reshaped so its thickness 
equals the minimum, with its concentration adjusted to conserve
ice volume.

Small amounts of either open water or ice area can be created due
to numerical diffusion associated with horizontal advection. To 
reduce the possibility of roundoff error corrupting the ice state
owing to very small amounts of sea ice, any ice
category whose area is less than an adjustable minimum (typically
$5 \times 10^{-6}$) is added to the nearest ice filled category,
if one is available. Small amounts of open water (typically less
than $1 \times 10^{-6}$) are eliminated by increasing ice
concentration equivalently in the thinnest ice category. Any
remaining small ice areas in a grid box are set to zero in such
a way that ice and snow volume are conserved across the
entire hemisphere's ice pack (i.e. renormalization factors are applied
across the entire hemispheric ice pack to compensate exactly for
setting small ice areas and volumes to zero).

\vskip 8pt
\noindent {\bf 4.2 Thermal Properties}

The ice area ($A_n$) and volume ($V_n$) were introduced in the Section 2.
The ice internal energy ($E_n$) is proportional to the ice volume: 
$$
E_n = q_n V_n   \eqno(14)
$$
where the proportionality function $q_n$ is termed the energy of melting,
or enthalpy. $q_n$
is the internal energy of the ice per unit volume. It is derived from the 
basic thermodynamic relation between the applied heat $Q$ for the given heat 
capacity of sea ice $c_i$ and the resulting temperature change from $T$ to $T^\prime$:
$$
Q = \int_T^{T^\prime} \rho_i c_i dT     \eqno(15)
$$
where the ice density $\rho_i$ is a constant (see List of Physical Constants). Treating 
ice density as a constant is a limitation of the model. During the melt season, a layer 
of deteriorated ice 5-10 cm thick is often observed at the top surface, with a density 
of 500 kg m$^{-3}$ or less. Beneath this deteriorated ice, multiyear ice contains 
air-filled pores that can reduce its density to 700-800 kg m$^{-3}$ in the upper 
30-50 cm (Bitz, 2000). A constant ice density implies that all drained brine 
pockets are filled with melt water and not air.

The storage of latent heat in brine pockets is accounted for explicitly by 
using the heat capacity of Bitz and Lipscomb (1999), originally from Ono (1967):
$$
c_i(T,S) = c_0 + {{L_i \mu S} \over {T^2}}     \eqno(16)
$$
where $c_0$ (J kg$^{-1}$ $^\circ$C$^{-1}$) is the heat capacity for pure water ice, 
$L_i$ (J kg$^{-1}$) is the latent heat of fusion of ice, $S$ (ppt) is the ice 
salinity, $T$ ($^\circ$C) is the temperature, and $\mu$ ($^\circ$C ppt$^{-1}$) 
is the empirical constant in the melting temperature ($T_{melt}$) and salinity 
relation:
$$ 
         T_{melt} = - \mu S.     \eqno(17)
$$
(See List of Physical Constants.)  For each category $n$, 
Eq. 15 is evaluated using
this heat capacity from temperature $T$ to the melting temperature $T_{melt}$
at salinity S:
$$
         q_n(T,S) = - \rho_i c_0 (T_{melt} - T) - \rho_i L_i 
         \left( 1 + {{\mu S} \over {T}} \right).     \eqno(18)
$$

The enthalpy $q_n$ is defined to be negative, so that $\vert q_n \vert$ 
is the amount of energy required to melt a unit volume of sea ice of salinity 
$S$ and temperature $T$. With this sign convention, a positive amount of heat 
$\vert q_n(T,S) \vert$ must be applied to raise the ice temperature from $T$ to 
$T_{melt}$ at salinity $S$, resulting in a rise of internal energy from $E_n<0$ 
to zero.

For snow, the heat required to change its temperature below melting is small compared 
to the latent heat of fusion, and thus for simplicity is ignored. 
Further, snow is fresh and therefore has zero salinity. Hence, the 
amount of energy required to melt a unit volume of snow is given by:
$$
         q_s = - \rho_s L_i      \eqno(19)
$$
where $\rho_s$ is the constant snow density. 
Note that since $q_s$ is a constant, the snow internal energy is proportional to $V_{sn}$,
so that an explicit snow internal energy is not a state variable. When required, the snow 
internal energy is computed from:
$$
E_{sn} = q_s V_{sn}.       \eqno(20)
$$
The snow/ice surface temperature is an important quantity for determining the heat and mass 
exchange between the atmosphere and the snow/ice surface. The surface temperature $T_{sn}$ 
for the $n^{th}$ category varies rapidly with changing forcing conditions, and because it
is used as an initial condition for the thermodynamic surface energy calculation, it
is treated as a state variable. The area-weighted surface temperature is used for 
conservation and transport, in both thickness and physical space.

Snow and ice thickness and ice temperature are not state variables, but they
can be diagnosed as follows. Ice and snow thickness are computed from the ice 
area and volume and snow volume, respectively, as:
$$
      h_n = V_n/A_n \quad h_{sn} = V_{sn}/A_n .      \eqno(21)
$$
Ice temperature can be diagnosed from the energy of melting. As discussed in Section 
4.6 on vertical heat conduction, the ice is divided vertically
into a number of layers. For each layer there is an internal energy and volume. From these
an energy of melting ($q_n$) can be computed for each layer (Eq. 14), and
hence a layer temperature from the solution to the quadratic equation (Eq. 18)
$$
         \rho_i c_0 T^2  
         - (q_n + \rho_i c_0 T_{melt} + \rho_i L_i)T - \rho_i L_i \mu S = 0.       \eqno(22)
$$
The solution yields one temperature below $T_{melt}$ and another above $T_{melt}$,
which is discarded.

\vskip 8pt
\noindent {\bf 4.3 Input from the Coupler}

Fields received by the ice model from the coupler are shown in 
Table 2. They include atmospheric/oceanic states 
and fluxes. Atmospheric states must be available to the ice model because
with an ice thickness distribution, it is necessary for the ice
model, rather than for the coupler (as in CSM1), to compute fluxes over each
ice thickness category and aggregate them. 

All of the fields received are on the T-grid (see Section 3.4).
However, the vector fields of surface wind, surface ocean current and tilt
are projected onto geographical latitude-longitude directions. These vectors
are first rotated to the displaced pole grid directions using a T-grid rotation 
angle($\chi^t$) calculated from a U-grid rotation angle ($\chi^u$) provided 
by a grid input dataset. The rotated surface wind is then used
on the T-grid to calculate atmosphere/ice fluxes, including stresses (see
Section 4.5). The resulting atmosphere/ice stress, as well as the ocean surface 
current and tilt, are then bilinearly interpolated with area weights to the
U-grid for use in the dynamics (see Section 4.10).

Specifically, let ($u_g^t,v_g^t$) represent a vector field of components ($u,v$),
where the subscript ``$g$'' refers to geographic, and the superscripts ``$t$''
and ``$u$'' to the T-grid and U-grid respectively. Similarly, let a subscript 
``$dp$'' refer to the displaced pole grid. All 
the vector fields received from the coupler are then ($u_g^t,v_g^t$) fields.
For these to be useful on the displaced pole grid, they must first be rotated
as follows:
$$
\eqalign{
      (u_{dp}^t)_{ij}  &= (u_g^t)_{ij} cos(\chi^t_{ij}) + (v_g^t)_{ij} sin(\chi^t_{ij}) \cr
      (v_{dp}^t)_{ij}  &= -(u_g^t)_{ij} sin(\chi^t_{ij}) + (v_g^t)_{ij} cos(\chi^t_{ij}) \cr
}                  \eqno(23)
$$
where ($ij$) are the longitude/latitude indices of the displaced pole grid. The 
atmosphere winds in the ($u_{dp}^t,v_{dp}^t$) form can be used to directly compute 
atmosphere/ice stresses. However, these stresses, as well as the ocean currents
and tilts, are required to be on the U-grid for the dynamic calculation. Therefore,
the following interpolation from the T-grid to the U-grid is required:
$$
\eqalign{
      (u_{dp}^u)_{ij}  = {1 \over 4}(A^t_{ij}(u_{dp}^t)_{ij} + A^t_{i+1j}(u_{dp}^t)_{i+1j} 
      + A^t_{ij+1}(u_{dp}^t)_{ij+1} + A^t_{i+1j+1}(u_{dp}^t)_{i+1j+1}) / A^u_{ij} \cr
      (v_{dp}^u)_{ij}  = {1 \over 4}(A^t_{ij}(v_{dp}^t)_{ij} + A^t_{i+1j}(v_{dp}^t)_{i+1j} 
      + A^t_{ij+1}(v_{dp}^t)_{ij+1} + A^t_{i+1j+1}(v_{dp}^t)_{i+1j+1}) / A^u_{ij} \cr
}      \eqno(24)
$$
where $A^t_{ij}$ is the T-grid box area, and $A^u_{ij}$ is the U-grid 
box area.

The freezing/melting potential $F_{Qoi}$ is calculated
in the ocean model and received by the ice model as input. In the ice model,
three forms of ice are distinguished: {\bf{frazil}} (which forms directly in
the ocean surface layer), {\bf{congelation}} (which forms at the ice base), and
{\bf{snow-ice}} (which forms by flooding of snow-covered ice). Frazil ice formation 
is determined by the ocean model, and the other two by the ice model.

If the ocean surface layer 
temperature ($T_o$) falls below freezing (at fixed temperature $T_{of}$), frazil
ice forms such that the heat flux $F_{Qoi}$ restores the ocean temperature
to freezing:
$$         
F_{Qoi} = \rho_o c_o h_o (T_{of}-T_o) / \Delta t       \eqno(25)
$$
where    $\rho_o c_o$  is the product of ocean density and heat capacity,
         $h_o$           is the surface layer thickness,
         $\Delta t$          is the coupling time step for the ocean (i.e. the time
between exchanges of data with the coupler, usually one day).
If ($T_o < T_{of}$) then $F_{Qoi} > 0$ and frazil ice forms (note that all CCSM fluxes 
are positive downwards).

For completeness, we discuss the salinity rejected in frazil ice formation, although
it has no effect on the sea ice model. The salinity adjustment $\Delta S$ in the 
ocean model due to brine rejection is:
$$    
\Delta S = (S_o - S_i) F_{Qoi} \Delta t / (\rho_o h_o L_i)       \eqno(26)
$$
where  $S_o$  is the constant reference ocean salinity,
       $S_i$  is the constant reference salinity of sea ice,
       $\rho_o$  is the constant ocean density, and
       $L_i$  is the latent heat of fusion of sea ice. 
Note that the ocean freezing temperature $T_{of}$ is kept constant independent 
of salinity, and the latent heat of fusion $L_i$ is 
identical in the ocean and sea ice models.

The total downwards shortwave flux is the sum of the four components, which
are defined in Table 2: 
$$    
F_{SWDN}=F_{SWvsdr}+F_{SWvsdf}+F_{SWnidr}+F_{SWnidf}.        \eqno(27)
$$

For the rest of the document, $\Delta t$ will represent both the physical and
the coupling time step for the sea ice model, which for CSIM5 is one hour.
It is possible to run the model with a physical time step for the dynamics
smaller than the coupling time step.

\vskip 8pt
\noindent {\bf 4.4 Snow and Ice Albedo}

Snow and ice albedos are important for computing the absorption of shortwave 
radiation in the snow/ice system, and hence snow/ice albedo feedback 
(Curry et al., 1995). The physics of this absorption and scattering is very 
complex (Ebert and Curry, 1993; Grenfell et al., 1994), but here it is simplied 
significantly. The snow and ice albedo formulas are basically those of CCSM2, 
which have been lowered somewhat to yield better sea ice simulation in CCSM3.
The albedo depends upon spectral band, snow thickness, ice thickness and 
surface temperature. 

Snow and ice spectral albedos (visible $=vs$, wavelength $< 0.7\mu m$ and 
near-infrared $=ni$, wavelength $> 0.7\mu m$) are distinguished, as both snow and 
ice spectral reflectivities are significantly higher in the $vs$ band than in the 
$ni$ band. This two-band separation represents the basic spectral dependence. Thus, 
we ignore the near-infrared spectral structure, with generally decreasing 
reflectivity with increasing wavelength (Ebert and Curry, 1993).

The zenith angle dependence of snow and ice is ignored (Ebert and Curry, 1993;
Grenfell et al., 1994), and therefore the distinction between downwelling direct 
(dr) and diffuse (df) shortwave radiation. The error in this approximation is 
probably no larger than .05 . Solar elevation angles in the polar regions are low.
The zenith angle dependence for the albedo affects the downwelling direct radiation,
while clear skies are relatively rare in polar regions. Horizontal variations in 
snow/ice topography are also ignored, which affect scattering and transmission into 
the surface through shadowing effects and through variations in the angle of the 
surface above the horizon. 

Snow albedo depends strongly on snow age (i.e. grain size, Grenfell et al., 1994), 
and on surface temperature (i.e. melting or non-melting conditions, (Ebert and Curry, 
1993). Sea ice albedo depends on ice thickness (Allison et al., 1993), as well as 
the presence of melt ponds (Ebert and Curry, 1993). In addition, snow only partially 
covers a surface if there are strong topographic variations (Allison et al., 1993).

Here we ignore the dependence of snow albedo on age, but retain the melting/non-melting 
distinction and thickness dependence. Dry snow spectral albedos for the $n^{th}$ 
category are: 
$$
\eqalign{
  \alpha_{vsdfn}^s(dry) = & {0.96}  \cr
  \alpha_{nidfn}^s(dry) = & {0.68}  \cr}   \eqno(28)
$$ (Note that these are less than the corresponding CCSM2 albedos at 0.98 and 0.70
respectively.) These values are roughly consistent with those of Grenfell et al. 
(1994) and Ebert and Curry (1993). In the case of the measurements of Grenfell et al. 
(1994), these dry snow albedos are slightly lower than clear sky values for low sun 
and for limited cloud cover, corresponding to the spring-time high values prior to 
significant melt (Curry et al., 2001). These albedos are only slightly higher than 
those for late summer conditions with early snow fall under cloudy skies.

To represent melting snow albedos, the surface temperature is used. Springtime 
warming produces a rapid transition from sub-zero
to melting temperatures, while late fall values transition more slowly to
sub-zero conditions. This is approximated by a temperature dependence out 
to $-1^\circ$C. Let $T_{snc}$ represent the snow/ice surface temperature
for category $n$ in $^\circ$C. If $T_{snc} = T_{sn}-T_{melt} \ge -1^\circ$C then
$$
\eqalign{
  \Delta T_s  =  & {T_{snc} + 1.0}    \cr
  \alpha_{vsdfn}^s(melt) = & {\alpha_{vsdfn}^s(dry) - 0.10 \Delta T_s}  \cr
   \alpha_{nidfn}^s(melt) = & {\alpha_{nidfn}^s(dry) - 0.15 \Delta T_s}  \cr}
                                          \eqno(29)
$$
The lowest albedos at $0^\circ$C are .86 and .53 for visible and near-ir
respectively, lower by about .02 than those of Ebert and Curry (1993) (the
corresponding CCSM2 albedos are .88 and .56 respectively). If the surface 
temperature $T_{snc} < -1^\circ$C, the dry snow albedos are used; otherwise 
the melt albedos.

For bare non-melting sea ice, albedo depends on thickness and spectral
band. If $h_n < 0.5$ m then
$$
\eqalign{
  \alpha_{vsdfn}(dry) = & \alpha_o(1-fh) + \alpha_{vsdfn}(thick)fh  \cr
  \alpha_{nidfn}(dry) = & \alpha_o(1-fh) + \alpha_{nidfn}(thick)fh  \cr }
                                          \eqno(30)
$$
where $\alpha_o$ is the open ocean diffuse albedo (as for snow and sea ice,
the zenith angle dependence of the ocean surface is also ignored),
$$
  fh = min( tan^{-1}(c_{fh}~h_i)/tan^{-1}(c_{fh}~0.5) , 1.0)        \eqno(31)
$$
$c_{fh}=4$ (for CCSM2, $c_{fh}=5$), and the thick, non-melting sea ice albedos are:
$$
\eqalign{
  \alpha_{vsdfn}(thick) = & {0.73}  \cr
  \alpha_{nidfn}(thick) = & {0.33}  \cr}         \eqno(32)
$$
which are the asymptotic values for ice thicker than $0.5$ m (corresponding
CCSM2 values are 0.78 and 0.36 respectively). These expressions represent a 
rough fit to the data of Allison et al. (1993), with the limiting cases for 
zero ice thickness that of the open ocean albedo $\alpha_o$, and for ice thicker 
than 0.5 m that of the thick ice case of Ebert and Curry (1993). The inverse 
tangent functional form approximates the theoretical dependence of ice albedo 
on thickness.

For bare melting sea ice, melt ponds can significantly lower the area 
averaged albedo. This effect is crudely approximated by the following
temperature dependence.
If $T_{snc} \ge -1^\circ$C, where $T_{snc} = T_{sn} - T_{melt}$, then
$$
\eqalign{
  \alpha_{vsdfn}(melt) = & {\alpha_{vsdfn}(dry) - 0.075 \Delta T_s }  \cr
  \alpha_{nidfn}(melt) = & {\alpha_{nidfn}(dry) - 0.075 \Delta T_s }  \cr}
                                                 \eqno(33)
$$
This results in minimum spectral albedos of .655 and .255 for visible
and near-ir respectively, or a rough broad band albedo (summertime spectral 
ratios of visible and near-ir of .53 and .47 respectively) of .467 (corresponding
spectral albedos for CCSM2 were .705 and .285 respectively, and a broadband
of .508). As for the case of snow, if the surface temperature $T_{snc} < -1^\circ$C, 
the dry sea ice albedos are used; otherwise the melt albedos.

The horizontal fraction of surface covered with snow is 
$$
       f_{sn} = {{h_{sn}} \over {h_{sn}} + 0.02}            \eqno(34)
$$ 
This expression is approximately in keeping with snow depth dependence of 
albedo from Ebert and Curry (1993), and from measurements of albedo on snow 
covered Antarctic sea ice Allison et al. (1993). We arrived at the value of 
$.02$ m in the denominator to achieve a good agreement with SHEBA data. This
means that $2$ cm of snow will cover $50\%$ of the horizontal sea-ice area;
the rest is assumed to be bare sea-ice.

Combining ice and snow albedos by averaging over the horizontal coverage
results in
$$
\eqalign{
  \alpha_{vsdfn} = & {\alpha_{vsdfn}(1 - f_{sn}) + f_{sn} \alpha_{vsdfn}^s}  \cr
  \alpha_{nidfn} = & {\alpha_{nidfn}(1 - f_{sn}) + f_{sn} \alpha_{nidfn}^s}  \cr}
                                                         \eqno(35)
$$
As noted above, the direct albedos are assumed identical to the diffuse.
These formulas are limited for thin, bare melting sea-ice to be greater 
than the ocean albedo $\alpha_o$.

This crude albedo, when compared with SHEBA measurements (Curry et al., 2001),
is able to approximately represent the major albedo regimes of 
springtime pre-melt dry snow, melting snow cover, dry bare ice, bare ice with 
melt ponds, and early fall freeze with light snow. However, the parameterization
is arbitrary and inconsistent across the snow/ice thermal and physical state, 
angle, spectral, and snow/ice thickness dependencies. It does not allow for
any consistent addition of organic and/or non-organic materials, and is not
consistent with internal heating and penetration into the underlying ocean.
There is no allowance for snow aging nor for explicit melt ponds. While a first
step, there is much room for improvement in this albedo parameterization. 

For diagnostic purposes, it is useful to have an aggregate broad band 
surface albedo for the history file (see Section 6):
$$
  \alpha_{bb} = .29\alpha_{vsdr} + .24\alpha_{vsdf} + .31\alpha_{nidr} + .16\alpha_{nidf}
                                              \eqno(36)
$$
The relative weights are only rough estimates of typical surface flux in each spectral
band and incident angle. A broad band albedo consistent with the downwelling shortwave
fluxes listed in Table 2 and the total shortwave flux in Eq. 27 could easily be computed
in future versions.

\vskip 8pt
\noindent {\bf 4.5 Ice to Atmosphere Flux Exchange}

Atmospheric states and downwelling fluxes, along with surface states and
properties, are used to compute atmosphere-ice shortwave and longwave fluxes, 
stress, sensible and latent heat fluxes. Surface states are temperature $T_{sn}$ and
albedos $\alpha_{vsdrn}$, $\alpha_{vsdfn}$, $\alpha_{nidrn}$, $\alpha_{nidfn}$
(see Section 4.4), while surface properties are longwave emissivity
$\varepsilon$ and aerodynamic roughness $z_i$ (note that these properties in 
general vary with ice thickness, but are here assumed constant). Additionally, 
certain flux temperature derivatives required for the ice temperature calculation 
are computed, as well as reference diagnostic surface air temperature and
specific humidity. 

The following formulae are used for the $n^{th}$ category for the absorbed 
shortwave fluxes and upwelling longwave flux:
$$       F_{SWvsn} = F_{SWvsdr}(1 - \alpha_{vsdrn}) + F_{SWvsdf}(1 - \alpha_{vsdfn}) \eqno(37) $$
$$       F_{SWnin} = F_{SWnidr}(1 - \alpha_{nidrn}) + F_{SWnidf}(1 - \alpha_{nidfn}) \eqno(38) $$
$$       F_{SWn} = F_{SWvsn} + F_{SWnin} \eqno(39) $$
$$       F_{LWUPn} = -\varepsilon\sigma_{sb} T_{sn}^4+(1-\varepsilon)F_{LWDN} \eqno(40) $$
The downwelling shortwave flux and albedos distinguish between visible ($vs, \lambda < 0.7\mu m$),
near-infrared ($ni, \lambda > 0.7\mu m$), direct ($dr$) and diffuse ($df$) radiation for each category. 
Note that the upwelling longwave flux has a reflected component from the downwelling longwave
whenever $\varepsilon < 1$.

For the $n^{th}$ category stress components, sensible and latent heat flux, the following bulk formulae 
are used (NCAR CSM Flux Coupler, 1996):
$$       \tau_{axn} = \rho_a r_{mn} u_n^\ast u_a   \eqno(41) $$
$$       \tau_{ayn} = \rho_a r_{mn} u_n^\ast v_a   \eqno(42) $$
$$       F_{SHn} = \rho_a c_a r_{hn} u_n^\ast \left( \theta_a - T_{sn} \right) \eqno(43) $$
$$       F_{LHn} = \rho_a L_s r_{en}  u_n^\ast \left( q_a - q_s (T_{sn}) \right) \eqno(44) $$
where:
$$       q_s (T_{sn}) = (q_1/ \rho_a) e^{-q_2/T_{sn}} \eqno(45) $$
$$       c_a = C_p ( 1 + C_{pvir} q_s(T_{sn})) \eqno(46) $$
$$       C_{pvir} = (C_{pwv}/C_p) - 1 . \eqno(47) $$
$L_s$ in Eq. 44 is the latent heat of sublimation from ice to vapor.
$q_s(T)$ is the surface saturation specific humidity for either ice or ocean at temperature $T$ 
in Kelvins (the values of $q_1,q_2$ for ice were kindly supplied by Xubin Zeng of
the University of Arizona), $C_p$ is the specific heat of dry air and $C_{pwv}$ of water vapor
(see List of Physical Constants for values of constants). The exchange coefficients for momentum, 
sensible and latent heat for each category are $r_{mn}$, $r_{hn}$, and $r_{en}$ respectively. 

The bulk formulae are based on Monin-Obukhov similarity theory. Among boundary layer scalings,
this is the most well tested (Large, 1998). It is based on the assumption that in the surface
layer (typically the lowest tenth of the atmospheric boundary layer), but away from the 
surface roughness elements, only the distance from the boundary and the surface kinematic fluxes
are important in the turbulent exchange. The fundamental turbulence scales that are formed from 
these quantities are the friction velocity $u_n^\ast$, the temperature and moisture fluctuations
$\theta_n^\ast$ and $q_n^\ast$ respectively, and the Monin-Obukhov length scale $L_n$:
$$       u_n^\ast      = r_{mn} V_{mag}    \eqno(48) $$
$$       \theta_n^\ast = r_{hn} (\theta_a - T_{sn})    \eqno(49) $$
$$       q_n^\ast      = r_{en} (q_a - q_s(T_{sn}))    \eqno(50) $$
$$       L_n           = u_n^{\ast 3} / (\kappa F_n)    \eqno(51) $$
with
$$       V_{mag} = max(1.0,\sqrt{u_a^2 + v_a^2}),  \eqno(52) $$
to prevent zero or small fluxes under quiescent wind conditions, $\kappa$ is 
von Karman's constant (0.4), and $F_n$ is the bouyancy flux, defined as:
$$       F_n = {{u_n^\ast} \over {g_e}}
\{ {{\theta_n^\ast} \over {\theta_{vn}}} + {{q_n^\ast} \over {z_v^{-1}+q_a}} \} \eqno(53) $$
with $g_e$ the gravitational acceleration and the virtual potential temperature $\theta_v 
= \theta_a(1+z_vq_a)$ where $z_v=\rho_{wv}/\rho_a - 1$.

Similarity theory holds that the vertical gradients of mean horizontal wind, potential
temperature and specific humidity are universal functions of stability parameter
$\zeta = z / L$, where $z$ is height above the surface ($\zeta$ is positive for a 
stable surface layer and negative for an unstable surface layer). These universal similarity
functions are determined from observations in the atmospheric boundary layer 
(Hogstrom, 1988) though no single form is widely accepted. Integrals of the 
vertical gradient relations result in the familiar logarithmic mean profiles, from 
which the exchange coefficients can be defined, where $\zeta_n = z_a / L_n$:
$$ r_{mn} = r_0 \{1+{{r_0} \over {\kappa}}(ln(z_a/z_{ref})-\chi_m(\zeta_n))\}^{-1}
\eqno(54) $$
$$ r_{hn} = r_0 \{1+{{r_0} \over {\kappa}}(ln(z_a/z_{ref})-\chi_h(\zeta_n))\}^{-1}
\eqno(55) $$
$$ r_{en} = r_{hn} \eqno(56) $$
with the neutral coefficient $r_0$ over ice:
$$ r_0 = {{\kappa} \over {ln(z_{ref}/z_i)}}, \eqno(57) $$
and over ocean:
$$ r_0 = (.0027/V_{mag} + .000142 + .0000764V_{mag})^{1/2} \eqno(58) $$
where $z_{ref}$ is a reference height (presently $10$ m). Note that the
square of $r_0$ is often referred to as the neutral drag coefficient $C_d$, and
that for the sea ice aerodynamic roughness value in the List of Physical Constants, 
$C_d = r_0^2 = 1.6 \times 10^{-3}$. The flux profile functions (integrals of the 
similarity functions mentioned above) for momentum m and heat/moisture h are:
$$   \chi_m(\zeta_n) = \chi_h(\zeta_n) = -5\zeta_n \eqno(59) $$
for stable conditions ($\zeta_n > 0$). For unstable conditions ($\zeta_n < 0$):
$$   \chi_m(\zeta_n) =  ln \{(1+X_n(2+X_n))(1+X_n^2)/8\} - 2 tan^{-1}(X_n) 
+ 0.5\pi  \eqno(60) $$
$$   \chi_h(\zeta_n) = 2 ln\{(1+X_n^2)/2\} \eqno(61) $$
with
$$   X_n = (max((1-16\zeta_n)^{1/2}),1.)^{1/2}. \eqno(62) $$

The stability parameter $\zeta_n$ is a function of the turbulent scales and 
thus the fluxes, so an iterative solution is necessary. The coefficients are
initialized with their neutral value $r_0$, from which the turbulent scales,
stability, and then flux profile functions can be evaluated. This order is
repeated for five interations to ensure convergence to an acceptible solution.

The surface temperature derivatives required by the ice temperature calculation are
evaluated as:
$$ {{\partial F_{LWUPn}} \over {\partial T_{sn}}} = -4\varepsilon\sigma_{sb} T_{sn}^3  \eqno(63) $$
$$ {{\partial F_{SHn}} \over {\partial T_{sn}}}   = - \rho_a c_a r_{hn} u_n^\ast \eqno(64) $$
$$ {{\partial F_{LHn}} \over {\partial T_{sn}}}   = - \rho_a L_s r_{en}  u_n^\ast 
{{\partial q_s(T_{sn})} \over {\partial T_{sn}}} \eqno(65) $$
where the small temperature dependencies of $c_a$, the exchange coefficients $r_{hn}$ and
$r_{en}$ and velocity scale $u_n^\ast$ are ignored.

For diagnostic purposes, an air temperature ($T_{REFn}$) and specific humidity
($Q_{REFn}$) are computed at the reference height of $z_{2m}=2m$. We define
$\zeta=\pm z_{2m}/z_a$, where the $\pm$ refers to stable/unstable conditions
respectively. Using Eq. 62 along with Eq. 59 (stable) and Eq. 61 (unstable), we
define $\chi_2$. Hence for temperature:
$$ f_{int} = (r_{hn} / \kappa)(-ln(z_{2m}/z_{ref})-\chi_{hn}+\chi_2) \eqno(66) $$
$$ T_{refn} = \theta_a + (\theta_a - T_{sn})f_{int} \eqno(67) $$
$$ T_{refn} = T_{refn} - .01z_{2m} \eqno(68) $$
For specific humidity we have:
$$ f_{int} = (r_{en} / \kappa)(-ln(z_{2m}/z_{ref})-\chi_{hn}+\chi_2) \eqno(69) $$
$$ Q_{refn} = q_a + (q_a - q_{sn})f_{int}. \eqno(70) $$

\vskip 8pt
\noindent {\bf 4.6 Vertical Heat Conduction}

Vertical heat conduction follows Maykut and Untersteiner (1971) and Bitz and
Lipscomb (1999). This section is also drawn from Bitz (2000) with modifications 
to match the notation used in this report, and with minor modifications for 
assumptions unique to this model.

To represent the vertical transfer of heat through the ice, we allow the ice internal energy 
(Eq. 14) to vary with level $z$, where $z$ is vertical depth measured positive 
downwards from the ice surface. The number of layers of ice ($L$) in each category
is four, with each layer thickness $\Delta h_n =h_n/L$ where $h_n$ is from Eq. 21.
The internal energy for each layer can be solved for an equivalent layer temperature
(Eq. 22). Vertical heat transfer is then calculated for $(l=1...L)$ vertical layers in the 
sea ice and one layer of overlying snow (if present). A staggered vertical grid 
is used, with temperature and salinity defined at layer midpoints and conductivity defined 
at layer interfaces. Layers at the top and bottom are referred to as surface layers, and 
those away from the surfaces as interior layers. In this section, the superscript is reserved 
for the time index $m$, and the category index $n$ is implied; the subscript $s$ on $T$ 
denotes the surface and the subscript $0$ denotes the snow layer. 

The invariant vertical salinity profile is represented by
$$
         S(w) = 1.6 \left\{ 1 - \cos \left( \pi w^{{{0.407} \over {0.573 + w}}} \right) \right\}
\eqno(71) $$
with the normalized coordinate $w$ calculated for each category as
$$ w = z/h_n,  0 \le w \le 1 \eqno(72) $$ Thus the salinity profile is independent of category. 
The profile varies from $0$ ppt at ice surface to $3.2$ ppt at ice base; salinity values 
at $w$ of 0.25, 0.50 and 0.75 are 1.5, 2.8 and 3.1 ppt respectively. The profile is meant 
to represent multi-year ice for which surface melt has reduced salinity near the top 
layers with respect to the lower. Note that this salinity profile need not be consistent
with the reference salinity used in ice-ocean salt exchange ($S_i$; see the List of 
Physical Constants). The vertically averaged salinity from the above profile is 2.3 ppt.

The heat content change over the time interval $t$ to $t'$ corresponding to temperatures 
$T$ and $T'$, respectively, allowing for temperature dependent heat capacity, thermal conduction 
and internal absorption of penetrating solar radiation, is given by (see Eqs. 15,18):
$$
  \int^{T'}_T \rho_i c_i dT = 
   \rho_i c_0(T'-T) \left(1 + {L_i\mu S \over  c_0 T' T} \right) 
  =\int^{t'}_t \left( { \partial  \over
  \partial z} k { \partial T \over \partial z} + Q_{SW} \right) dt
\eqno(73) $$
where $c_i$ is from Eq. 16, $Q_{SW}$ is the absorbed shortwave flux, 
and the thermal conductivity $k$ is either that for snow or ice. For snow, $k = k_s$ 
is a constant, while for ice:
$$
    k(S,T) = k_{fi} + {{\beta S} \over {T}}
\eqno(74) $$
where $k_{fi}$ and $\beta$ are empirical constants from Untersteiner (1961).
$Q_{SW}$ is given by:
$$
  Q_{SW} = - {{d} \over {dz}} \{ I_{0vs} e^{-\kappa_{vs} z} + I_{0ni} e^{-\kappa_{ni} z} \}
\eqno(75) $$
and $I_{0vs},I_{0ni}$ are the fractions of absorbed shortwave radiation in the visible 
and near-ir that penetrate the surface, respectively, given by (see Eq. 37):
$$
I_{0vs} = 0.70 F_{SWvsn} (1-f_{sn})
\eqno(76) $$
$$
I_{0ni} = 0.0
\eqno(77) $$
where $f_{sn}$, the horizontal fraction of surface covered by snow, is given by 
Eq. 34. It is assumed that no shortwave radiation penetrates the snow covered 
surface. The ice spectral extinction coefficients $\kappa_{vs}$ and $\kappa_{ni}$ 
are for the visible and near-infrared bands respectively (Gary Maykut, personal 
communication). To compute the penetration factors (.70 and .0) for the visible and 
near-ir radiation respectively, a surface layer of $~5$ cm thick was assumed. However, 
for the surface energy balance calculation (see Section 4.6.1, Eq. 83) the surface 
layer thickness is not explicitly used. Note that there is no distinction made 
between direct or diffuse shortwave: in effect, we assume shortwave radiation 
penetrating the surface is diffuse. Also, the functional form of penetration in
Eq. 75 is that for pure absorption; in effect, we assume all scattering occurs
at the surface while the rest of the ice is purely absorbing.

The heat equation (Eq. 73) is discretized using a
backwards-Euler, space-centered scheme. Using the staggered grid with
$T_l$ representing the layer temperature and $k_l$ representing
conductivity at the layer interfaces, for interior layers we have
$$ 
 \rho_i c_0(T_l^{m+1}  -  T_l^m)
  \left(1 + {L_i\mu S_l \over  c_0 T_l^{m+1} T_l^m} \right) = 
   { \Delta t  \over  \Delta h^m  } \left( 
    k_{l+1}^m { T_{l+1}^{m+1} - T_l^{m+1}  \over  \Delta h^m  } 
    - k_l^m   { T_l^{m+1} - T_{l-1}^{m+1}  \over  \Delta h^m  } 
    +  I^m_l \right),   
\eqno(78) $$
where $\Delta h^m =h^m/L$, the conductivity is
$$
k_l^m = k\left({ S_l+S_{l+1} \over 2 } ,
         { T_l^m+T_{l+1}^m \over 2 }  \right),
\eqno(79) $$
and the absorbed solar radiation is
$$
I_l^m  = I_{0vs} (e^{-\kappa_{vs} l \Delta h^m} - 
                  e^{-\kappa_{vs} (l+1) \Delta h^m}) +
         I_{0ni} (e^{-\kappa_{ni} l \Delta h^m} - 
                  e^{-\kappa_{ni} (l+1) \Delta h^m}).
\eqno(80) $$

For a purely implicit backward scheme, $k$ should be evaluated at the
$m+1$ time level. However, when $k$ is evaluated at time level $m$,
experiments show that the solution is stable and converges to the same
solution one gets when evaluating $k$ at $m+1$.

The discrete heat equation for the surface layers is modified slightly
from Eq. 78 to maintain second-order accuracy for $\partial
T/\partial z$. The equation for the bottom layer ($l=L$) is
$$
\eqalign{ 
 \rho_i c_0&(T_L^{m+1}- T_L^m)
  \left(1 + {L_i\mu S_L \over  c_0 T_L^{m+1} T_L^m} \right) =\cr 
  & { \Delta t  \over  \Delta h^m  } 
    \left(
    3 k_{L+1}        { T_b - T_L^{m+1} \over \Delta h^m    } -
 {1 \over 3} k_{L+1} { T_b - T_{L-1}^{m+1} \over \Delta h^m}
   - k_L^m             { T_L^{m+1} - T_{L-1}^{m+1}  \over  \Delta h^m  } 
          +  I^m_L \right), \cr}
\eqno(81) $$
where the $L+1$ interface in contact with the underlying ocean is assumed to
be at temperature $T_b=T_{of}$ and salinity $S_b$ (given by Eqs. 71,72 with
$w=1$), and therefore conductivity $k_{L+1} = k(S_b,T_b)$.
The equations for the top surface depend on the surface conditions,
of which there are four possibilities, as outlined in Table 5. If the 
snow depth is too small, numerical solutions are unreliable, and hence 
the insulating effects of snow are ignored for depths below $h_{smin}$.

\vskip 15pt
\vbox{
\hang \noindent {{\it Table 5.} 
Top Surface Boundary Cases}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 & Snow Accumulated & Melting \cr
\noalign{\hrule height 1pt}
      case I    &   yes        &       no      \cr
      case II   &    no        &       no      \cr
      case III  &   yes        &       yes     \cr
      case IV   &    no        &       yes     \cr
\noalign{\hrule height 1pt}
}}}
\vskip 10pt

\vskip 8pt
\noindent {\bf 4.6.1 Surface Boundary Conditions}

Water fluxes from the atmosphere are received by the ice model in the 
form of rain $F_{RN}$ and snow $F_{SNW}$ (see Table 2). Presently, the rain 
is assumed to run off into the ocean without modification to the snow and 
ice, and adds to the water flux into the ocean $F_{Wo}$ (see Table 3). 
Unless otherwise noted, all fluxes have an implied time step index $m$.

Snow (see Eq. 21) is assumed to accumulate on the surface as:
$$
h_{sn}^{m+1/2} = h_{sn}^m + F_{SNW}\Delta t / \rho_s .  \eqno(82)   
$$

The boundary condition at the ice surface results from a balance of fluxes
$$
   F_{TOPn}(T_{sn}) = F_{SWvsn}-I_{0vs} + F_{SWnin}-I_{0ni}  
   + \varepsilon F_{LWDN} - \varepsilon 
   \sigma_{sb} T_{sn}^{4} + 
   F_{SHn} + F_{LHn} + k_n {{dT_n} \over {dz}} {\bigg\vert}_{z=0}, \eqno(83)   
$$
where $z$ is the vertical depth measured downwards from the ice surface.

If $F_{TOPn} (T_{sn}) \le 0$ then no surface melting occurs. If 
$F_{TOPn} (T_{sn}) > 0$ then $T_{sn} = T_{melt}$ and surface melting 
when snow is present proceeds according to
$$
        F_{TOPn}(T_{sn}) = q_s {{dh_{sn}} \over {dt}}.    \eqno(84)   
$$
When no snow is present surface melting proceeds according to
$$
        F_{TOPn}(T_{sn}) = q_n(T_n,S) {{dh_n} \over {dt}},    \eqno(85)   
$$
where $T_n$ and $S$ are the ice temperature and salinity of the top ice
layer, respectively. The boundary condition at the ice bottom is
$$
      - F_{BOT} - k_n {{\partial{T_n}} \over {\partial{z}}} = q_n(T_n,S) 
        {{dh_n} \over {dt}}    \eqno(86)   
$$
where $F_{BOT}$ is the heat flux from the ocean (see section 4.8, Eq. 131),
$k_n$ is the ice thermal conductivity, and $T_n$,$S$ are the temperature and
salinity of the bottom ice layer, respectively. 

Four cases must be distinguished to solve for vertical heat conduction, 
based on presence or absence of snow, and whether melting or non-melting
conditions pertain.

\vskip 8pt
\noindent {\bf 4.6.2 Case I: Snow accumulation with no melting}

The discrete heat equation for the uppermost snow layer (recall that 
the subcript $0$ on T denotes the snow layer) is: 
$$ 
 \rho_s c_s(T_0^{m+1} - T_0^m) 
  =  { \Delta t  \over  h_s^m }  
\left[ k_1^m { T_1^{m+1} - T_0^{m+1} \over (\Delta h^m+h_s^m)/2} 
   - \alpha k_s { T_0^{m+1} - T_s^{m+1} \over h_s^m  } 
      - \beta k_s { T_1^{m+1} - T_s^{m+1} \over h_s^m  } 
         \right].       \eqno(87)    
$$
The heat equation solver is formulated for the general case where the heat
capacity of snow $c_s$ may be specified, although here is set to $0$.
The parameters $\alpha$ and $\beta$ are defined to give second-order accurate
spatial differencing for $\partial T / \partial z$ across the
changing layer spacing at the snow/ice boundary;
$$
\eqalign{ 
  \alpha & = { h_s^m + \Delta h^m/2 \over h_s^m/2 }~{2 \over h_s^m + \Delta h^m} h_s^m\cr
  \beta  & = {-h_s^m/2 \over h_s^m + \Delta h^m/2 }~{2 \over h_s^m + \Delta h^m} h_s^m.\cr 
}     \eqno(88)    
$$
The conductivity at the snow--ice interface is found by equating
conductive fluxes above and below the interface;
$$
k_1^m={2 k_s k(S_1,T_1^m) \over h_s^m k(S_1,T_1^m) + \Delta h^m k_s} ~
   {h_s^m+\Delta h^m \over 2}.     \eqno(89)    
$$
Because $T_s$ is below melting, a flux boundary condition is used, and an additional 
equation is required in the coupled set:
$$ 
F_o(T_s^{m+1})+\alpha k_s { T_0^{m+1} - T_s^{m+1} \over h_s^m }
      +\beta k_s { T_1^{m+1} - T_s^{m+1} \over h_s^m }= 0,     \eqno(90)    
$$
where $F_o(T_s^{m+1})$ is the sum of all terms on the right-hand side
of Eq. 83 except $k \partial T / \partial z$.  The net surface
flux $F_o(T_s^{m+1})$ is approximated as linear in $T_s^{m+1}$; thus
$$
F_o(T_s^{m+1})\sim F_o(T_s^m)+\left.{\partial F_o \over \partial
    T_s}\right|_{T_s^m}(T_s^{m+1}-T_s^m).     \eqno(91)    
$$
with
$$
\left.{\partial F_o \over \partial T_s}\right|_{T_s^m} =
\left.{\partial F_{LWUP} \over \partial T_s}\right|_{T_s^m} +
\left.{\partial F_{SH} \over \partial T_s}\right|_{T_s^m} +
\left.{\partial F_{LH} \over \partial T_s}\right|_{T_s^m}     \eqno(92)    
$$
(see Eqs. 63-65).

To simplify our set of equations, we define
$$
\hat c_l^{m+1}= \rho_i \left( c_0 + {L_i \mu S \over T_l^{m+1} T_l^m}\right),     \eqno(93)    
$$
where the hat implies that $\hat c_l^{m+1}$
depends on $T_l^m$ as well as on $T_l^{m+1}$,
and
$$
\chi_l^{m+1}= {\Delta t \over \Delta h^m } {1 \over \hat c_l^{m+1}}.     \eqno(94)    
$$
Also, let
$$
k_l= {k_l^m \over \Delta h^m}.     \eqno(95)    
$$
for $l \ge 2$ and
$$
k_0= {k_s \over h_s^m}     \eqno(96)    
$$
$$
k_1= {k_1^m \over (\Delta h^m + h_s^m)/2}     \eqno(97)    
$$
and suppress the index $m$ for $I_l^m$, so that 
for interior layers ($l=1...L-1$),
$$
T_l^{m+1}-T_l^m=\chi_l^{m+1}\left[k_{l+1}(T_{l+1}^{m+1}-T_l^{m+1})
               -k_l(T_l^{m+1}-T_{l-1}^{m+1})+I_l\right]     \eqno(98)    
$$
and at the bottom layer
$$
\eqalign{
       T_L^{m+1}-T_L^m&=\chi_L^{m+1}\left[3k_b(T_b-T_L^{m+1})
                           -{1 \over 3}k_b(T_b-T_{L-1}^{m+1}) \right. \cr
                 & \ \ \ \ \ \ \ \ \ \ \ \ \left.   
                            -k_L(T_L^{m+1}-T_{L-1}^{m+1})+I_L\right] \cr
}     \eqno(99)    
$$
where $k_b = k_{L+1} / \Delta h^m$.
The equation describing the snow layer is written
$$
\rho_s c_s (T_0^{m+1}-T_0^m)={\Delta t \over h_s^m } 
                  \left[ k_1(T_1^{m+1}-T_0^{m+1})
                  -\alpha k_0(T_0^{m+1}-T_s^{m+1})
                  -\beta  k_0(T_1^{m+1}-T_s^{m+1}) \right].     \eqno(100)    
$$
%where $\rho_s c_s$ remains in the numerator to allow the option of
%setting $c_s=0$. 
Finally, the flux boundary condition becomes
$$
      F_o(T_s^m)+ \left. {\partial F_o \over \partial T_s} \right|_{T_s^m}
      (T_s^{m+1}-T_s^m) =
 -\alpha k_0(T_0^{m+1}-T_s^{m+1}) -\beta k_0(T_1^{m+1}-T_s^{m+1}).  \eqno(101)    
$$

The complete set of coupled equations for case I can be written with
all of the terms that explicitly depend on temperature at the $m+1$
time step gathered on the right-hand side:

$$
\eqalign{
      -F_o(T_s^m)+ \left.{\partial F_o \over 
               \partial T_s}\right|_{T_s^m}T_s^m &=
           T_s^{m+1}\left(\left.{\partial F_o \over \partial T_s}
             \right|_{T_s^m}-\alpha k_0 - \beta k_0 \right) \cr
         & \hskip 2.0truepc + T_0^{m+1} \alpha k_0 + T_1^{m+1} \beta k_0 \cr
      \rho_s c_s T_0^m  &=
            T_s^{m+1}\left( -{\Delta t \over h_s^m } \right) (\alpha k_0 
                  +  \beta k_0) \cr
         & \hskip 2.0truepc  + T_0^{m+1}\left( \rho_s c_s 
          + {\Delta t \over h_s^m }  (\alpha k_0+k_1) \right) \cr
         & \hskip 2.0truepc  + T_1^{m+1}{\Delta t \over h_s^m } (\beta k_0-k_1) \cr
      T_l^m + \chi_l^{m+1} I_l &=
            T_{l-1}^{m+1}(-\chi_l^{m+1} k_l) \cr
         & \hskip 2.0truepc + T_l^{m+1}(1+\chi_l^{m+1} k_l + \chi_l^{m+1} k_{l+1} ) \cr
         & \hskip 2.0truepc + T_{l+1}^{m+1}(-\chi_l^{m+1} k_{l+1})  \cr
      T_L^m + \chi_L^{m+1} I_L + {8 \over 3} \chi_L^{m+1} k_bT_b &=
            T_{L-1}^{m+1}\left(-{1 \over 3}\chi_L^{m+1} k_b -
            \chi_L^{m+1} k_L \right) \cr
         & \hskip 2.0truepc + T_L^{m+1}(1+3 \chi_L^{m+1} k_b + \chi_L^{m+1} k_L ). \cr
}     \eqno(102)    
$$
These equations are subsequently related to 
the following abbreviated form
$$
\eqalign{
      r_s &= T_s^{m+1} b_s + T_0^{m+1}c_s+ T_1^{m+1}d_s \cr
      r_0 &= T_s^{m+1} a_0 + T_0^{m+1}b_0+ T_1^{m+1}c_0 \cr
      r_1 &= T_0^{m+1} a_1 + T_1^{m+1}b_1+ T_2^{m+1}c_1 \cr
            & \vdots \cr
      r_L &= T_{L-1}^{m+1} a_L + T_L^{m+1}b_L. \cr
}     \eqno(103)    
$$
The first two rows can be combined to eliminate the coefficient on
$T_1^{m+1}$ in the first row, allowing the set to be written in
tridiagonal form:
$$
 \matrix{
{ r = \left[ \matrix{
                r_s c_0 - r_0 d_s \cr
                r_0 \cr
                r_1 \cr
               \vdots     \cr 
                    }  \right] }
&\hskip 1.0truepc
 A= \left[ \matrix{
                 b_s c_0 - a_0 d_s & c_s c_0 - b_0 d_s &      &     \cr
                 a_0               &  b_0              & c_0  &     \cr
                                   &  a_1              & b_1  & c_1 \cr
                                   &                   &      & \ddots  \cr
                  } \right]
&\hskip 1.0truepc
 T = \left[ \matrix{
               T_s^{m+1} \cr
               T_0^{m+1} \cr
               T_1^{m+1} \cr
               \vdots     \cr
                   } \right].
}     \eqno(104)    
$$
Because the matrix A depends on $\chi_l^{m+1}$, which in turn depends
on $T_l^{m+1}$, the system of equations is solved iteratively. An initial
guess is used for the temperature dependence of $\chi_l^{m+1}$,
and then $\chi_l^{m+1}$ is updated successively after each iteration.
Under most conditions the method approaches a solution in less than
four iterations with a maximum error tolerance of $\Delta T_{err}$ for
$T_l$ with an initial guess of $T_l^{m+1}=T_l^m$.
% end of case I

\vskip 8pt
\noindent {\bf 4.6.3 Case II: Snow free with no melting}

Nearly the same method applies when the ice is snow free, except one
less equation is needed to describe the evolution of the temperature
profile. The equation for the uppermost ice layer is written
$$ 
\eqalign{
     \rho_i c_0(& T_1^{m+1} - T_1^m)  
  \left(1 + {L_i\mu S_1 \over  c_0 T_1^{m+1} T_1^m} \right) \cr 
  &= { \Delta t  \over  \Delta h^m  } 
   \left( 
      k_2^m { T_2^{m+1} - T_1^{m+1}  \over \Delta h^m } 
      -3 k_1^m { T_1^{m+1} - T_s^{m+1} \over \Delta h^m  } 
    +{1 \over 3} k_1^m { T_2^{m+1} - T_s^{m+1} \over \Delta h^m } 
          +  I^m_1 \right),  \cr}     \eqno(105)    
$$
where $k_1^m=k(S_1,T_1^m)$. After the definitions from Eqs. 93--95
are applied, Eq. 105 becomes
$$
T_1^{m+1}-T_1^m=\chi_1^{m+1}\left[k_2(T_2^{m+1}-T_1^{m+1})
                             -3k_1(T_1^{m+1}-T_s^{m+1})
                     +{1 \over 3}k_1(T_2^{m+1}-T_s^{m+1})+I_1^m\right].    \eqno(106)    
$$
The flux boundary condition follows
after linearizing $F_o(T_s^{m+1})$ in $T_s^{m+1}$:
$$
     F_o(T_s^m)+\left.{\partial F_o \over \partial T_s}\right|_{T_s^m}
     (T_s^{m+1}-T_s^m) =
    -3k_1(T_1^{m+1}-T_s^{m+1}) +{1 \over 3} k_1(T_2^{m+1}-T_s^{m+1}).    \eqno(107)    
$$
The complete set of coupled equation includes Eqs. 102 for
layers 2 to L with the following two equations for the surface and
upper ice layer:
$$
\eqalign{
      -F_o(T_s^m)+\left.{\partial F_o \over \partial T_s}\right|_{T_s^m} T_s^m &=
            T_s^{m+1}\left(\left.{\partial F_o \over \partial T_s}\right|_{T_s^m}
          -k_1 {8 \over 3}\right) + T_1^{m+1}3 k_1 + T_2^{m+1}( -k_1/3 ) \cr
      T_1^m + \chi_1^{m+1} I_1^m &=
            T_s^{m+1}\left(-\chi_1^{m+1} k_1 {8 \over 3}\right) \cr
          & + T_1^{m+1}(1+\chi_1^{m+1} k_2 + 3 \chi_1^{m+1} k_1 ) \cr
          & + T_2^{m+1}(-\chi_1^{m+1} k_2 - {1 \over 3}\chi_1^{m+1} k_1), \cr
}    \eqno(108)    
$$
which can be written
$$
\eqalign{
      r_s &= T_s^{m+1} b_s + T_1^{m+1}c_s+ T_2^{m+1}d_s \cr
      r_1 &= T_s^{m+1} a_1 + T_1^{m+1}b_1+ T_2^{m+1}c_1. \cr
}    \eqno(109)    
$$
These two equations can be combined to eliminate the coefficient on
$T_2^{m+1}$, allowing the set to be written in tridiagonal form:
$$
\matrix{
 r = \left[ \matrix{
                r_s c_1 - r_1 d_s \cr
                r_1 \cr
                r_2 \cr
               \vdots     \cr
                   }     \right]
&\hskip 1.0truepc
 A= \left[ \matrix{
                 b_s c_1 - a_1 d_s & c_s c_1 - b_1 d_s &      &     \cr
                 a_1               &  b_1              & c_1  &     \cr
                                   &  a_2              & b_2  & c_2 \cr
                                   &                   &      & \ddots  \cr
                  } \right]
&\hskip 1.0truepc
 T = \left[ \matrix{
               T_s^{m+1} \cr
               T_1^{m+1} \cr
               T_2^{m+1} \cr
               \vdots     \cr
                   }     \right].
}    \eqno(110)    
$$
As for case I, this system of equations must be solved iteratively. 
% end of II

\vskip 8pt
\noindent {\bf 4.6.4 Case III: Snow accumulation with melting}

Case III describes melting conditions in the
presence of a snow layer at the surface.  Here a temperature boundary
condition is used, which simplifies the solution because the first row 
in Eqs. 102 is not needed and $T_s=T_{melt}$ in the second row.  
Hence the complete set of coupled equations is identical to Eqs. 102
for layers 1 to L, with the addition of an equation for the snow layer,
$$
   \rho_s c_s T_0^m + T_{melt} {\Delta t \over h_s } 
      (\alpha + \beta)k_0 =
            T_0^{m+1}\left[ \rho_s c_s + {\Delta t \over h_s }
                   (k_1+\alpha k_0 ) \right]
          - T_1^{m+1} {\Delta t \over h_s } (k_1-\beta k_0).  \eqno(111)
$$
This set of equations can be written in tridiagonal form, without the
need to eliminate any terms, as was required in cases I and II. 
However, the solution must still be iterated.

\vskip 8pt
\noindent {\bf 4.6.5 Case IV: No snow with melting}

Like case III, case IV describes melting conditions, but here the sea 
ice is snow free. Hence, the first two rows of Eqs. 102
are not needed, and $T_s=T_{melt}$ for $l=1$.  The
set of coupled equations comprises those from Eqs. 102 for
layers 2 to L and the following equation for layer 1:
$$
      T_1^m + \chi_1^{m+1} I_1^m + T_{melt} \chi_1^{m+1} k_1 {8 \over 3} = 
          T_1^{m+1}\left(1+\chi_1^{m+1} k_2 + 3 \chi_1^{m+1} k_1 \right)
        + T_2^{m+1}\left(-\chi_1^{m+1} k_2 - 
          {1 \over 3}\chi_1^{m+1} k_1 \right).  \eqno(112)
$$
As in case III, this set of equations can immediately be written in the
tridiagonal form and solved iteratively.

\vskip 8pt
\noindent {\bf 4.7 Thickness Changes: Top, Bottom and Snow-to-Ice Conversion}

The energy of melting of snow is a constant, given by Eq. 19, 
while the energy of melting for each layer $l$ of sea ice depends on the 
temperature and salinity ($T_l$ and $S_l$, respectively) of the layer 
according to Eqs. 17 and 18:
$$
q_l= - \rho_i c_0(-\mu S_l - T_l) - \rho_i L_i \left( 1 +{\mu S_l \over T_l} \right).
     \eqno(113)    
$$
The energy balance at the top and bottom surfaces determines the top 
melt and bottom melt and growth rates of the sea ice. From the top surface 
flux balance in Eq. 83, if $F_{\rm TOP}(T_{sn}) > 0$, then the upper surface 
is fixed at the melting temperature and $F_{\rm TOP}$ is used for melting,
according to
$$
\left. \delta h_s \right|_{\rm melt} = { F_{\rm TOP}(T_{melt}) \Delta t
    \over q_s }     \eqno(114)    
$$
where $ \left. \delta h_s \right|_{\rm melt}$ is the change in snow thickness
due to melt. For bare ice:
$$ 
\left. \delta h \right|_{\rm melt} = { F_{\rm TOP}(T_{melt}) \Delta t \over q_1 } ,
     \eqno(115)    
$$
where $ \left. \delta h \right|_{\rm melt} $ is the change in the top layer
thickness due to ice melt, and $q_s$ and $q_1$ are the energy of melting of 
the snow and the top layer of the ice, respectively. If all the snow is melted
in one time step, any residual heat is applied to melt the ice. Snow and ice melt 
water is assumed to drain to the ocean below without any effect on the intervening 
snow and ice.

Sublimation of snow occurs when $F_{LH}<0$ (regardless of $T_s$), according to
$$ 
  \left. \delta h_s \right|_{\rm sub} = 
  - { F_{LH} \Delta t  \over (q_s-\rho_s L_v) }     \eqno(116)    
$$
where $ \left. \delta h_s \right|_{\rm sub} $ is the change in 
snow thickness due to sublimation, and for ice according to
$$ 
  \left. \delta h\right|_{\rm sub} = 
  - { F_{LH} \Delta t \over (q_1-\rho_i L_v) }     \eqno(117)    
$$
if the ice is snow free, where $ \left. \delta h\right|_{\rm sub} $ is the change
in the top layer ice thickness due to sublimation. The same set of equations applies
when $F_{LH}>0$ for condensation on the ice or snow.

The bottom-surface energy balance is (see section 4.8)
$$ 
  \left. \delta h\right|_{\rm basal} = 
     { (- F_{BOT} - k {\partial T \over \partial z}) \Delta t \over
      q_b },     \eqno(118)    
$$
where $ \left. \delta h\right|_{\rm basal} $ is the change in the lowest ice
layer thickness due to basal freezing or melting, $F_{BOT}$ is the heat supplied
to the ice from the underlying ocean (see section 4.8, Eq. 131), the basal
heat conduction is 
$$ 
k {\partial T \over \partial z} =
    3 k_b        { T_b - T_L^{m+1} \over \Delta h^m    } -
 {1 \over 3} k_b { T_b - T_{L-1}^{m+1} \over \Delta h^m}     \eqno(119)    
$$
accurate to second order, where the subscript ``b'' refers to ice base
(i.e. $k_b = k_{L+1}$; see discussion following Eq. 81), and
$$
q_b=\cases{ q_L; & $\left. \delta h \right|_{basal}<0$ \cr
 - \rho_i c_0(-\mu S_b - T_b) - \rho_i L_i \left( 1 +{\mu S_b \over
  T_b} \right); & $\left. \delta h \right|_{basal}>0$.}     \eqno(120)    
$$
$\left. \delta h\right|_{\rm basal} > 0$ represents formation of
congelation ice, and is treated as a (negative) water flux to
the ocean; while $\left. \delta h\right|_{\rm basal} < 0$ represents
basal ice melting which is added to the water flux to the ocean
(see section 4.13).

Snow to ice conversion is allowed. This occurs if the
snow layer overlying the sea ice becomes thick enough to
depress the snow-ice interface below freeboard (the ocean
surface). The interface height is:
$$ 
z_{int} = h - (\rho_s h_s + \rho_i h) / \rho_o.       \eqno(121)    
$$
If $z_{int} < 0$, then an amount of snow in depth equal to $-z_{int} \rho_i / \rho_s$
is removed from the snow layer and added to the ice. It is 
assumed that ocean water floods the depressed snow, and then
converts it into ice of thickness $-z_{int}$. The energy of melting
of the newly formed ice is: $q_{flood} = q_s \rho_i / \rho_s$. 

The energy of melting of the ice and snow layers needs to be adjusted
when the layer spacing changes after growth/melt, sublimation/condensation,
and flooding.  The adjusted energy of melting is
$$
q_l'= \cases { \sum_{k=1}^L w_{k,1} q_k - q_{flood} {z_{int} \over \Delta h'} ; & $l=1$ \cr
      \sum_{k=1}^L w_{k,l} q_k; & $1<l<L$ \cr
      \sum_{k=1}^L w_{k,L} q_k + q_b ~
     {\rm max}({\left. \delta h \right|_{\rm basal} \over \Delta
       h'},0); & $l=L$. }   \eqno(122)    
$$
where $w_{k,l}$ are weights computed from the relative overlap of 
layer $l$ with each layer $k$ from the old layer spacing and $\Delta
h'$ is the new layer spacing.  

As discussed in section 4.3, and further in section 4.13, for the purposes of
ocean-ice salt exchange, sea ice is assumed to have a constant reference salinity 
$S_i$. This implies that any ice formation or melt must be accompanied by a salt 
exchange with the ocean in order to maintain the constant ice reference salinity
(apart from frazil ice formation whose salinity exchange with the ocean is accounted 
for in the ocean model itself- see section 4.3). 

Let $\delta V = \delta h A + h \delta A$ be the change in ice volume due either
to a $\delta h$ change in ice thickness or a $\delta A$ change in ice area. The
following equations give the salt flux exchange with the ocean for ice melt,
net congelation (or basal) at ice base, net sublimation and condensation,
and snow ice formation respectively:
$$  
F_S(ice~melt) = +\rho_i S_i \delta V_{ice~melt} / \Delta t \eqno(123)   
$$  
$$  
F_S(cong~ice) = -\rho_i S_i \delta V_{cong~ice} / \Delta t \eqno(124)   
$$  
$$  
F_S(sub~cond) = -\rho_i S_i \delta V_{sub~cond} / \Delta t \eqno(125)   
$$  
$$  
F_S(snow~ice) = -\rho_i S_i \delta V_{snow-ice} / \Delta t \eqno(126)   
$$  
When ice melts, salinity is given to the ocean (positive flux), but for 
non-frazil ice formation salinity must be drawn out of the ocean. This 
is because the associated water flux exchange for non-frazil ice formation 
with the ocean is fresh.

\noindent {\bf 4.8 Lateral Formation, Side and Bottom Melt Fluxes}

Lateral formation and melt occurs depending on the sign of the
freezing/melting potential $F_{Qoi}$ (see Section 4.3, Eq. 25).

For {\bf{positive}} $F_{Qoi}$, frazil ice formation occurs.  
A volume of frazil ice $V_f = F_{Qoi}\Delta t/q_f$ is formed, where 
$q_f=-\rho_i L_i$ is the enthalpy of frazil ice formation.
If there is open water (i.e. $A_0>0$), then the thickness of the newly formed frazil 
ice is $h_f=V_f/A_0$. If $h_f < h_{min}$, new ice 
of thickness $h_{min}$ is assumed to form over area $A_0^\prime = V_f / h_{min}$, else 
$A_0^\prime = A_0$, where primes denote updated quantities. This new ice is added to 
category 1 so long as $h_f<h_1^\star$, the thickest ice allowed for that category. 
If the new ice is thicker than this limit, or if there is no open water (i.e. 
$A_0=0$), then the new ice is distributed evenly over all categories. The lower 
flux boundary condition is $F_{BOT}$=0.

For new frazil ice evenly distributed over all categories, 
$V_{surp} =~  {V_f - A_0 h_1^\star}$ is the surplus ice volume 
over category 1, after which $V_f =~ {A_0 h_1^\star}$, so that:
$$
\eqalign{
  V_n'       =~   & {V_n + V_{surp}}A_n/A \cr
  E_n'(z)    =~   & E_n(z) + (V_{surp}A_n/AL) q_f \cr}     \eqno(127)     
$$
where $A$ is the total ice concentration over all categories and $E_n$ 
is adjusted for each level.

For new frazil ice added to category 1, for which even distribution
of ice to thicker categories may have occured, we have:
$$
\eqalign{
  A_1'    =~   & {A_1 + A_0'} \cr
  V_1'    =~   & {V_1 + V_f } \cr
  T_{s1}' =~   & {(T_{s1} A_1 + T_{of} A_0')/A_1'} \cr
  E_1'(z) =~   & E_1(z) + (V_f / L) q_f \cr}     \eqno(128)     
$$
where $E_1'$ is adjusted for each level.

For {\bf{negative}} $F_{Qoi}$, heat is available to melt ice.  This flux
is partitioned into heat available for side melt and bottom melt based on first
assuming $F_{Qoi}$ is dominated by shortwave radiation, and then assuming shortwave
radiation absorbed in the ocean surface layer above the mean ice thickness causes
side melting and below it causes basal melting. For the mean ice thickness 
$\bar{h} = V/A$, where $V$ and $A$ are the total ice volume and area respectively 
(see Section 4.13):
$$f_{bot} = {Re^{-\bar{h}/\zeta_1} + (1-R)e^{-\bar{h}/\zeta_2}} \eqno(129) $$
$$f_{sid} = {1-f_{bot} }  \eqno(130) $$    
where $R = 0.68$, $\zeta_1=1.2~m^{-1}$, $\zeta_2=28~m^{-1}$ (Paulson and Simpson, 1977)
and $f_{bot}$ and $f_{sid}$ are the fractions of bottom and side melt flux available, 
respectively. Thus the maximum fluxes available for melt are $f_{bot} F_{Qoi}$ and 
$f_{sid} F_{Qoi}$. The actual amount used for bottom melting, $F_{BOT}$, is based on 
boundary layer theory (Mcphee, 1992).
$$ 
F_{BOT} = max(-\rho_o c_o c_h u^* \Delta T, f_{bot} F_{Qoi})     \eqno(131)     
$$
where the empirical drag coefficient $c_h$=0.006, and
$$ \Delta T = max(T_o - T_{of}, 0)     \eqno(132)     $$
$$ u^* = min\{ \sqrt{(\tau_{ox}^2 + \tau_{oy}^2)/\rho_o}, u_{min} \}   \eqno(133)     $$
with $u_{min}$ the minimum allowed skin friction velocity $u^*$, and
the ocean-ice stresses $\tau_{ox}$ and $\tau_{oy}$ are from Eq. 157.

The heat flux for lateral melt is the product of the vertically averaged, aggregate
energy of melting of snow and ice, $E_{tot}$, with the interfacial 
melting rate $M_a$ and the total floe perimeter $p_f$ per unit floe area $A_f$.
$$E_{tot} = \rho_s L_s V_s + \sum_{n=1}^N \sum_{n=1}^{L_n} q_{nl} V_n / L_n \eqno(134)$$ 
where $q_{nl}$ is the energy of melting for layer 
$l$, and $L_n$ is the number of layers of sea ice for category $n$.
The interfacial melting rate is taken from the empirical expression of Maykut and
Perovich (1987)
based on Marginal Ice Zone Experiment observations: $M_a = m_1 \Delta T^{m_2}$, where
$m_1=1.6\times10^{-6} $m s$^{-1}$ deg$^{m_2}$ and $m_2=1.36$. The lead-ice perimeter
depends on the ice floe distribution and geometry. For a mean floe diameter $d$ and
number of floes $n_f$, $p_f = n_f \pi d$ and the floe area $A_f = n_f \eta_{lm} d^2$
(Rothrock and Thordike, 1984; Bitz 2000). As the heat flux for lateral melt is
$E_{tot}(p_f/A_f)M_a$, the actual amount used is:
$$ 
F_{SID} = max({{E_{tot} \pi} \over {\eta_{lm} d}} m_1 \Delta T^{m_2}, 
f_{sid} F_{Qoi})     \eqno(135)     
$$
where $\eta_{lm}=0.66$ (Rothrock and Thordike, 1984). Based partially on tuning 
and partially on the results of floe distribution measurements, the mean floe 
diameter of $d$=300 m was chosen. The ice area, volume, snow volume, and ice 
energy are all reduced by side melt in time $\Delta t$ by the fraction 
$R_{side} = \vert {{F_{SID}\Delta t} \over {E_{tot}}} \vert$.

The heat flux available that is actually used by the ice model is:
$$ 
F_{Qio} = F_{BOT} + F_{SID}    \eqno(136)     
$$

\vskip 8pt
\noindent {\bf 4.9 Linear Remapping}

In this section we evaluate the first right-hand side term in the distribution
equation (Eq. 1) due to thermodynamic processes, which can be thought
of as a transport in thickness space:
$$ 
         {{\partial g} \over {\partial{t}}} = 
       - {{\partial} \over {\partial h}} ({\dot{h}}g)      \eqno(137)   
$$ 
where ${\dot{h}}=dh/dt$. 
To evaluate this term we use the linear remapping method of Lipscomb (2001).
The method is similar to the 1D version of the incremental 
remapping algorithm of Dukowicz and Baumgardner (2000). 

The linear remapping method uses the integral form of the distribution equation.
Integrating the above transport equation between $h_{n-1}^\ast(t)$ 
and $h_n^\ast(t)$, where the boundary thicknesses $h_n^\ast,\{n=0,1,2...N\}$ 
are time dependent (i.e. following the motion), results in
$$ 
      \int_{h_{n-1}^\ast (t)}^{h_n^\ast (t)}   {{\partial g} \over {\partial{t}}} dh + 
\int_{h_{n-1}^\ast (t)}^{h_n^\ast (t)} {{\partial} \over {\partial h}} ({\dot{h}}g) dh = 0 .
     \eqno(138)   
$$ 
Evaluating the second term yields
$$ 
      \int_{h_{n-1}^\ast (t)}^{h_n^\ast (t)} {{\partial} \over {\partial h}} ({\dot{h}}g) dh = 
g(h_n^\ast) \left.{{dh^\ast} \over {dt}}\right|_n - g(h_{n-1}^\ast) \left.{{dh^\ast} \over {dt}}
\right|_{n-1} .     \eqno(139)   
$$ 
The integrated transport equation can be rewritten as
$$ 
         {{d} \over {dt}} \int_{h_{n-1}^\ast (t)}^{h_n^\ast (t)} g(h) dh = 0 ,  \eqno(140)   
$$ 
or using Eq. 2
$$ 
         {{d A_n} \over {dt}} = 0 .     \eqno(141)   
$$ 
This equation can be interpreted in a Lagrangian sense as a conservation
equation, where the time dependent limits are the boundaries of the Lagrangian
volume following the motion in thickness space. Differentiating the volume
($V_n = A_n h_n$) with time we have
$$ 
         {{d V_n} \over {dt}} = A_n {{dh_n} \over {dt}} .     \eqno(142)   
$$ 

The linear remapping method calculates the new Lagrangian boundaries, linearly 
approximates the distribution function, and transfers ice area and volume
to restore the original thickness boundaries. $g(h)$ is approximated by a 
series of linear piecewise continuous functions, and ice is transferred in 
small increments between categories.

The original ice thickness in category $n$ is $h_n^m$. After the thermodynamic
(both vertical and lateral) changes to ice thickness are computed, the new ice
thickness is $h_n^{m+1}$. Let the growth rate of ice thickness in category $n$ 
be $f_n = {\dot{h}_n} = (h_n^{m+1}-h_n^m)/ \Delta t$. The $m+1$ growth rates at 
$h_n^{\ast m}$ are estimated by interpolating between neighboring values of $f_n$:
$$
         f_n^{\ast} = f_n + {{(f_{n+1}-f_n)} \over
                             {(h_{n+1}^{m}-h_n^m)}}
                             (h_n^{\ast m}-h_n^{m})     \eqno(143)   
$$ 
The new boundary locations are then: 
$h_n^{\ast m+1} = h_n^{\ast m} + f_n^{\ast} \Delta t$. Note that in principle
the boundaries can shift any distance, but we require here that 
$h_{n-1}^{m+1} < h_n^{\ast m+1} < h_{n+1}^{m+1}$. 

If any category has no ice (i.e. $A_n=0$), while an adjacent category does
($A_{n+1} \not= 0$), the boundary is adjusted by the
same amount as the thickness in the non-zero category:
$h_n^{\ast m+1} = h_n^{\ast m} + f_{n+1}^{\ast} \Delta t$.

We approximate each Lagrangian volume linearly as $g_n(h) = g_{0n} + g_{1n} h$.
To evaluate the coefficients $\{g_{0n},g_{1n}\}$, we make use of the area
and volume constraints, the requirement of a positive distribution function, and the
minimum and maximum boundary conditions.
The area and volume constraints are
$$ 
         \int_{h_{n-1}^{\ast m} }^{h_n^{\ast m} } g(h) dh =
         \int_{h_{n-1}^{\ast m+1} }^{h_n^{\ast m+1} } g(h) dh = A_n^m     \eqno(144)   
$$ 
$$ 
         \int_{h_{n-1}^{\ast m+1} }^{h_n^{\ast m+1} } hg(h) dh = A_n^m h_n^{m+1}     \eqno(145)   
$$ 
Note that the area is conserved following the motion (i.e $A_n^m$ is a constant) 
but the volume changes to $A_n^m h_n^{m+1}$.

Let us redefine the $n^{th}$ category boundaries at time step $m+1$ as
$h_L=h_{n-1}^{\ast m+1}$ and $h_R=h_n^{\ast m+1}$ (henceforth
we drop the time and category indices). A positive linear function is constructed 
such that the area and volume integrals over the Lagrangian
volume are satisfied. We transform variables from $h$ to $\eta$
(a relative coordinate): for each Lagrangian volume $\eta = h - h_L$
and $\eta_n = h_n - h_L$ so that $g(\eta) = g_0 + g_1 \eta $, 
where $\eta$ ranges from $0$ to $\eta_R = h_R - h_L$. Note that 
the $h_n$ in $\eta_n=h_n-h_L$ is the $m+1$ value $h_n^{m+1}$.
$$
         A_n = \int_0^{\eta_R} g(\eta) d\eta = 
{{1} \over {2}}{\eta_R}^2 g_1 + \eta_R g_0     \eqno(146)   
$$ 
$$
         A_n h_n = \int_0^{\eta_R} (\eta+h_L) (g_1\eta + g_0) d\eta = 
h_LA_n + {{1} \over {3}} {\eta_R}^3 g_1 + {{1} \over {2}} {\eta_R}^2 g_0  \eqno(147)   
$$ 
We have two linear equations for $g_{0n}$ and $g_{1n}$ for the $n^{th}$
category
$$
{{1} \over {2}}{\eta_R}^2 g_{1n} + \eta_R g_{0n} = A_n  \eqno(148)   
$$ 
$$
{{1} \over {3}} {\eta_R}^3 g_{1n} + {{1} \over {2}} {\eta_R}^2 g_{0n} = A_n \eta_n \eqno(149)   $$ 
which have the solution
$$
g_{1n} = {{12A_n} \over {{\eta_R}^3}}(\eta_n - {{\eta_R} \over {2}})   \eqno(150)   
$$ 
$$
g_{0n} = {{6A_n} \over {{\eta_R}^2}}({{2\eta_R} \over {3}} - \eta_n).  \eqno(151)   
$$ 
Note that the sign of $g_{1n}$ is determined by $\eta_n - {{\eta_R} \over {2}} =
h_n - h_L - {{h_R-h_L} \over {2}} = h_n - {{1} \over {2}}(h_L+h_R)$. When $h_n$ is
greater than the Lagrangian midpoint, the slope is positive; when it is
less, it is negative.

As $g$ is linear, its maximum and minimum lie at the boundaries 
$\eta=0$ and $\eta=\eta_R$
$$
g(0) = {{6A_n} \over {{\eta_R}^2}}({{2\eta_R} \over {3}} - \eta_n) \eqno(152)
$$
and
$$
g(\eta_R) = {{6A_n} \over {{\eta_R}^2}}(\eta_n-{{\eta_R} \over {3}}) \eqno(153) .
$$
For $g(\eta)$ to be positive, both boundary values must be positive.
$g(0)$ is less than zero when $({{2\eta_R} \over {3}} - \eta_n) < 0$, or
$\eta_n > {{2\eta_R} \over {3}}$, and $g(\eta_R)$ is less than zero when 
$\eta_n < {{\eta_R} \over {3}}$, i.e. whenever $\eta_n$ lies outside the
central third of the Lagrangian thickness range. As just noted, whenever
$h_n$ is greater than the range mid point, the slope $g_{1n} > 0$; if it
is greater than $h_L + {{2} \over {3}}(h_R-h_L)$, the slope is so great
that the minimum value at $g(0)$ falls below zero; and conversely for
$h_n < h_L + {{1} \over {3}}(h_R-h_L)$, for which the slope becomes so
negative as to require $g(\eta_R)$ to be less than zero. 

For the case when $h_n$ falls in the first third of the Lagrangian thickness
range, we redefine the upper limit to $h_R^\prime$ as:
$h_n = h_L + {{1} \over {3}}(h_R^\prime-h_L)$ or
$h_R^\prime = 3h_n - 2h_L$ and set $g=0$ between $h_R^\prime$ and $h_R$.
Similarly, when $h_n$ lies in the upper third of the Lagrangian thickness 
range, the lower limit is redefined to $h_L^\prime$ as:
$h_n = h_L^\prime + {{2} \over {3}}(h_R-h_L^\prime)$ or
$h_L^\prime = 3h_n - 2h_R$ and set $g=0$ between $h_L$ and $h_L^\prime$.
In either case, the solutions for $g_{0n}$ and $g_{1n}$ can still be used as
long as the appropriate boundaries are defined.

Once $g(h)$ is constructed for each category, the thickness distribution is
remapped to the original category boundaries by transferring the appropriate
area $\Delta A_n$ and volume $\Delta V_n$. If the displaced boundary 
$h_n^{\ast m+1}$ has moved to the right, then:
$$
\eqalign{
     \Delta A_n^{m+1} = \int_{h_n^{\ast m}}^{h_n^{\ast m+1} } g_{n}(h) dh \cr
     \Delta V_n^{m+1} = \int_{h_n^{\ast m}}^{h_n^{\ast m+1} } hg_{n}(h) dh \cr} \eqno(154)
$$
If the displaced boundary has moved to the left, the limits of integration
are reversed.

The thinnest (i.e. $n=1$) and thickest (i.e. $n=N$) categories have special 
minimum and maximum boundary conditions respectively. For category 1, if
sea ice is growing in open water at a positive rate $f_0$, shift $h_1^\ast$
to the right by $f_0 \Delta t$. If sea ice is not growing in open water,
approximate the growth rate as $f_0 = f_1$. If $f_0<0$, reduce the ice area
by the integral of $g(h)$ from $0$ to $-f_0 \Delta t$, leaving the ice volume
fixed, as ice volume cannot cross the left boundary. For the right boundary,
$h_N^\ast$ varies with $h_N$. As $g(h)$ is linear, setting 
$h_N^\ast = 3h_N - 2h_{N-1}^\ast$ ensures $g(h_N^\ast) = 0$.

Snow volume, internal energy of ice and surface temperature are affected by thickness 
space transport as follows. Assuming that within each category snow depth varies linearly
with ice thickness, the snow volume transport is proportional to the ice volume 
transport (see Section 2.2).
The internal energy of ice is $E_n = q_n V_n$, where $q_n$ is the energy 
of melting of ice, a function of ice temperature and salinity. The internal 
energy of ice is proportional to the volume of ice, so that the new internal energy 
at $m+1$ is $E_n^{m+1}=q_n^{m+1}V_n^{m+1}$, where $q_n^{m+1}$ is the $m+1$ energy of 
melting after the vertical thermodynamic heat transfer has been computed. The surface 
temperature $T_{sn}^m$ changes with area due to thickness space transport.

\vskip 8pt
\noindent {\bf 4.10 Velocity}

Pack ice is composed of rigid ice floes, ranging in size from order $1$ m 
to greater than $10$ km. The characteristics of motion for the
pack ice are discontinuous slippage near shore, near rigid motion
under considerable wind forcing (i.e. nearly rate-independent stress),
small or zero tensile strength for both uniaxial and two dimensional
dilation, and high compressive strength.

To model this material, the resolved sea-ice is considered to be a 
highly fractured, closely packed, isotropic medium in which inter-floe
forces are contact stresses. The resolved spatial scales are considered
to be much larger than the scale of inhomogeneities (i.e. floes). We consider
only aggregate ice motion (across the ITD) in each grid box. Thus, we model sea-ice 
as a two dimensional continuum, whose momentum conservation is described by
$$
        {\bar m} {{\partial u} \over {\partial t}} = 
       -{\bar m} f k\times u + 
           \tau_a + \tau_o + 
          {\bar m}g_e\nabla{H_o} + \nabla\cdot\sigma     \eqno(155)    
$$ 
where $\bar m$ is the total mass of snow and ice per unit area given by
$$
          {\bar m} = \rho_s {\sum_{n=1}^N}V_{sn} + \rho_i {\sum_{n=1}^N}V_n ,  \eqno(156)    
$$ 
the non-linear $\bf{u}$ advection terms are ignored as they are negligibly 
small when the equations are scaled, $f$ is the Coriolis parameter, 
$\bf{k}$ is the local vertical unit vector, $\tau_a$ and $\tau_o$ are forces 
due to air and water stresses respectively, $g_e$ is the gravitational 
acceleration, $H_o$ is the sea surface slope and $\nabla\cdot\sigma$ is 
the force due to internal ice stress.

The air-ice stress $\tau_a$ is described in Section 4.5 (Eqs. 41,42 and following)
with aggregate $\tau_a = (1/A)\Sigma_{n=1}^N \tau_{an} A_n$. The ocean-ice stress 
$\tau_o$ is in the form of a nonlinear drag law
$$
           \tau_o = \rho_o c_d |\bf{u}_o - \bf{u}|(\bf{u}_o - \bf{u})     \eqno(157)    
$$ 
where $c_d$ is the ocean-ice drag coefficient, $\rho_o$ is the density of seawater and $\bf
{u}_o$  is the surface ocean current.  Note that while the drag coefficient $c_d$ varies 
with ice thickness for actual ice floes, here we assume it to be constant.

We note the following modification to the usual two dimensional continuum equation
used above for momentum conservation (Eq. 155), that better approximates point-wise
free drift theory. That theory states the actual (rather than cell-averaged) ice velocity
should be uniform under uniform forcing, regardless of ice area. As written, Eq. 155
will give different solutions under free drift (i.e. $\nabla\cdot\sigma = 0$) for
different sea ice areas but constant ice thickness (i.e. ${\bar m}$ is proportional to
sea ice area for constant thickness). The modification is to multiply the air-ice
and ocean-ice stresses by the concentration A, specifically, use ${\sum_{n=1}^N}A_n \tau_{an} 
+ A \tau_o$ in place of the $\tau_a + \tau_o$ terms in Eq. 155, so that the solutions to
the equation better approximate velocity for the sea ice itself, rather than cell-average.
For a complete description of this modification, see Hunke and Dukowicz (2003) and
Connolley et.al (2004).

The internal stress tensor $\sigma$ is a linear vector function of a vector 
argument, which gives the internal force on the material for a specified direction. 
The general constitutive law for sea-ice relates the stress to the rate of strain, and
generally includes elastic (linear and reversible) and plastic (non-linear and irreversible)
components. 

The CSM1 sea ice model assumes a cavitating fluid rheology, in which both elastic, 
shear and tensile stresses are ignored (Bettge et al. 1996), and the ocean-ice stress 
is linearized. The model suffers numerical grid convergence difficulties near the North 
Pole (Weatherly et al. 1998). The cavitating fluid rheology is useful in some circumstances, 
but is limited especially due to lack of shear stresses. A more realistic and generally 
accepted rheology is the viscous-plastic, or VP (Hibler 1979; Kreyscher et al., 2000).

The VP rheology derives from the general stress-strain relation for viscous fluids
$$
\sigma_{ij} = 2\eta {\dot\epsilon_{ij}} + (\zeta - \eta){\dot\epsilon_{kk}\delta_{ij}}
-{P \over 2}\delta_{ij}    \eqno(158)    
$$
where the stress tensor is $\sigma_{ij}$ ($i,j$ = component indices), 
the compressive strength $P$, $\delta_{ij}$ is the Kronecker delta
and the rate of strain tensor is
$$
           {\dot\epsilon_{ij}} = {{1} \over {2}} \left( {{{\partial{u_i}} \over {\partial{x_j}}} + 
                                 {{\partial{u_j}} \over {\partial{x_i}}}} \right),    \eqno(159)    
$$
the total linear rate of strain (divergence) is $\dot\epsilon_{kk} = \dot\epsilon_{11} +
\dot\epsilon_{22}$ and $\dot\epsilon_{12}=\dot\epsilon_{21}$ is the shear rate of strain
component. $\zeta$ and $\eta$ are bulk (i.e. linear) and shear viscosities respectively.
The general form of this plastic rheology satisfies the condition that the deformational
part of the plastic stress tensor vanishes for constant $\bf{u}$ and that for solid
body rotation no stress is produced. For the viscous flow to be dissipative both $\zeta$
and $\eta$ must be positive.

The {\bf{plastic}} assumption is that the flow obeys an idealized plastic behavior, 
namely, that stressed ice is motionless until a yield stress is obtained, after which 
the flow is irreversible and rate independent. The principal stress states for plastic 
deformation lie on the yield curve specified by a normalized convex yield function, 
while the (irreversible) deformation is given by the normal flow rule.

The VP rheology assumes an elliptical yield curve of specified ratio of major
to minor axes $e$. In terms of the stress and strain rate (or deformation) 
invariants, obtained by diagonalizing and transforming the symmetric stress 
($\sigma_I$,$\sigma_{II}$) and rate of strain tensors 
($\dot\epsilon_I$,$\dot\epsilon_{II}$) into pure compression (I)
and shear (II) components, we have the yield function 
$$ 
Y(\sigma_I,\sigma_{II}) = { { \left( \sigma_I + {{P} \over {2}}\right) ^2} \over
{\left( {{P} \over {2}}\right)^2}} +
                     {{\sigma_{II}^2} \over {\left( {{P} \over {2e}}\right)^2}} = 1 
    \eqno(160)    $$
which is chosen to lie in the lower left quadrant of principal stress space, 
corresponding to very small tensile stress but with finite compressional and shear 
stresses (the latter if $e$ is relatively small). With the normal flow rule
$\dot\epsilon_I = \lambda {\partial Y \over \partial \sigma_I}$,  
and $\dot\epsilon_{II} = \lambda {\partial Y \over \partial \sigma_{II}}$, 
the unknown $\lambda$ can be computed, resulting in
$$ \sigma_I =  \zeta \dot{\epsilon}_I - P/2     \eqno(161)    $$
$$ \sigma_{II} =  \eta \dot{\epsilon}_{II}     \eqno(162)    $$
$$ \zeta = P/2 \Delta     \eqno(163)    $$
$$ \eta =  \zeta/e^2     \eqno(164)    $$
where
$$
\dot\epsilon_I = {\partial u \over \partial x}+{\partial v \over \partial y}    \eqno(165)    
$$
$$
\dot\epsilon_{II} = \{
({\partial u \over \partial x}-{\partial v \over \partial y})^2 +
({\partial u \over \partial y}+{\partial v \over \partial x})^2
                    \}^{1 \over 2}     \eqno(166)    
$$
(Stern et al., 1995), and
$$
          \Delta = \left[ \dot{\epsilon}_I^2 + \dot{\epsilon}_{II}^2 / e^2 \right]^{1/2}.
    \eqno(167)    
$$
In the limit of zero strain rate (i.e. rigid solid with no deformation), the viscosities
become infinite. To regularize this behavior, the viscosities are bounded for sufficiently
small strain rates so that the sea ice moves as a linear viscous fluid undergoing slow 
creep. Minimum viscosities are set to prevent non-linear instabilities.

The {\bf{elastic viscous plastic}} (EVP) rheology (Hunke and Dukowicz 1997)
derives from the simple stress-strain 
relation for small strains: $\sigma_{ij} = E \epsilon_{ij}$, where $E$ (analogous to 
Young's modulus) is related to ice strength such that it increases as ice strength increases. 
Writing this relation in terms of rate of strain, and ignoring the non-linear advection as 
in the momentum equation, gives ${1 \over E} {\partial \sigma_{ij} \over \partial t} = 
\dot\epsilon_{ij}$. The stress tensor equation (Eq. 158) for viscous 
fluids can be solved in terms of the rate of strain, as
$$
\dot\epsilon_{ij} = {1 \over 2\eta} \sigma_{ij} + {(\eta - \zeta) \over 4\eta\zeta}
{\sigma_{kk}\delta_{ij}} + {P \over 4\zeta}\delta_{ij}
  \eqno(168)    $$
Combining these two rates of strain (elastic and plastic) to yield the total
rate of strain results in
$$
              {\dot\epsilon_{ij}} =
  {1 \over E}{{{\partial{\sigma_{ij}}}} \over {\partial{t}}} + 
  {1 \over 2 \eta} \sigma_{ij} + 
  { {\eta - \zeta} \over {4\eta\zeta}} \sigma_{kk} \delta_{ij} + 
  {P \over {4\zeta}} \delta_{ij}  \eqno(169)    
$$
In the limit $E \to \infty$ this rate of strain equation asymptotes to the pure VP rheology,
while for $\eta,\zeta \to \infty$, the purely elastic rheology is recovered. Hence, as
$\eta,\zeta \to \infty$ under conditions of very small strain rate, the elastic term controls
the solution behavior, and represents a regularization of the VP rheology.
Note that the elastic term in the stress tensor equation requires that the
stress tensor components become prognostic variables in EVP. This is in contrast to VP 
for which the stress tensor components are diagnostic.
The elastic parameter $E$ is given in terms of the bulk viscosity and a damping time scale
for elastic waves $T_{ew}$ and time step $\Delta t$ 
$$ E =  \zeta/T_{ew}   \eqno(170)    $$
$$ T_{ew} =  E_0 \Delta t.  \eqno(171)    $$
where $E_0$ is a constant less than 1.
The momentum and stress tensor equations are
$$
      {\bar m} { {\partial u} \over {\partial t}} + c^\prime u - {\bar m}fv =
           c^\prime u_o + \tau_{ax} -
           {\bar m}g_e{{\partial H_o} \over {\partial x}} + F_x
  \eqno(172)    $$
$$
      {\bar m} {{\partial v} \over {\partial t}} + c^\prime v + {\bar m}fu =
           c^\prime v_o + \tau_{ay} -
           {\bar m}g_e{{\partial H_o} \over {\partial y}} + F_y
  \eqno(173)    $$
$$
      { {\partial \sigma_{ij}} \over {\partial{t}} } + 
      { {e^2} \over {2T_{ew}}} \sigma_{ij} + 
      { {1 - e^2} \over {4T_{ew}}} \sigma_{kk} \delta_{ij} = 
      { {P} \over {2T_{ew}} \Delta^\prime} {\dot\epsilon_{ij}} - 
      { {P} \over {4T_{ew}}} \delta_{ij}
  \eqno(174)    $$
where $c^\prime = \rho_o c_d |\bf{u}_o - \bf{u}|$ and 
$\Delta^\prime = max(\Delta,\Delta_{min})$ and if $\Delta < \Delta_{min}$,
then $P$ is replaced by $P \Delta / \Delta_{min}$, which prevents residual ice 
motion due to spatial variations in $P$ for extremely small or zero rates of 
strain, where $\Delta_{min}=10^{-13} A^t$, $A^t$ is the T-grid area, section 3.4.
The stress divergence terms $F_x$ and $F_y$ are evaluated
for a general orthogonal curvilinear coordinate system subject
to the contraints that the discretization be dissipative
and includes grid curvature effects (see Hunke and Dukowicz, 2002).

It is convenient to introduce the divergence $D_D$, the horizontal tension
$D_T$ and shearing $D_S$ strain rates defined by:
$$ D_D = \dot{\epsilon}_{11} + \dot{\epsilon}_{22}   \eqno(175) $$
$$ D_T = \dot{\epsilon}_{11} - \dot{\epsilon}_{22}   \eqno(176) $$
$$ D_S = 2\dot{\epsilon}_{12}                        \eqno(177) $$
Letting $\sigma_1=\sigma_{11}+\sigma_{22}$ and $\sigma_2=\sigma_{11}-\sigma_{22}$, 
Eqs. 174 can be alternatively expressed as:
$$
{1\over E}{\partial\sigma_1\over\partial t} + {\sigma_1\over 2\zeta} 
  + {P\over 2\zeta} = D_D \eqno(178) $$
$$
{1\over E}{\partial\sigma_2\over\partial t} + {\sigma_2\over 2\eta} = D_T \eqno(179) $$
$$
{1\over E}{\partial\sigma_{12}\over\partial t} + {\sigma_{12}\over 
       2\eta} = {1\over 2}D_S .   \eqno(180)   $$
where we note that $\sigma_I = \sigma_1 / 2$, $\sigma_{II} = \sqrt{\sigma_2^2 / 4 + 
\sigma_{12}^2}$, $\dot{\epsilon}_I = D_D$, $\dot{\epsilon}_{II} = \sqrt{D_T^2 + D_S^2}$, 
the above definitions of $\zeta$ and $\eta$ are unchanged, while $\Delta$ is:
$$
\Delta = \left[D_D^2 + {1\over e^2}\left(D_T^2 + D_S^2\right)\right]^{1/2}.  \eqno(181) $$
Multiplying the momentum equations by $u$ and $v$ respectively and summing,
one can form the kinetic energy equation. Using the product rule,
$u_i {{\partial \sigma_{ij}} \over {\partial x_j}} = 
{{\partial} \over {\partial x_j}} (u_i \sigma_{ij}) -
\sigma_{ij} {{\partial u_i} \over {\partial x_j}} = 
{{\partial} \over {\partial x_j}} (u_i \sigma_{ij}) -
\sigma_{ij} \dot\epsilon_{ij}$,
where repeated indices imply summation over $(1,2)$. The area integral of the
${{\partial} \over {\partial x_j}} (u_i \sigma_{ij})$ term will vanish on lateral 
boundaries where $u_i=0$, or for open ocean where $\sigma_{ij}=0$. Thus, area 
integrals of the kinetic energy equation over the ice will result in dissipative 
internal stresses so long as D, defined by:
$$
D = \int\left(\sigma_{11}\dot\epsilon_{11} +
  2\sigma_{12}\dot\epsilon_{12} + \sigma_{22}\dot\epsilon_{22}\right) dA,  \eqno(182)
$$
is positive definite. Using the definitions of $D_D$, $D_T$ and $D_S$ we have:
$$
D = \int\left[{1\over 2}\sigma_1 D_D + {1\over 2}\sigma_2 D_T 
    + \sigma_{12}D_S \right]dA   $$
In steady state (${\partial \sigma \over \partial t} \rightarrow 0$), 
Eqs 178, 179 and 180 reduce to the viscous-plastic constitutive law. Using the 
steady state forms for $\sigma_1$, $\sigma_2$ and $\sigma_{12}$, $D$ can be 
written as:
$$
D = {1\over 2}\int P\left(\Delta-D_D\right)dA \ge 0,  \eqno(183) $$
which ensures that the work done by the internal stress is dissipative.
Note that $D$ is a scalar invariant independent of coordinate system.

To include grid curvature effects, consider an orthogonal, curvilinear 
coordinate system, and let $u$ and $v$ represent
velocity components along nondimensional coordinates $\xi_1$ and $\xi_2$
($0$ to $1$ across a grid box), with scale factors (the physical lengths 
of the grid box sides) $h_1$ and $h_2$, respectively. The rate of strain 
components are then:
$$
\dot{\epsilon}_{11} 
= {1\over h_1}\left({\partial u\over\partial\xi_1}
   +{v\over h_2}{\partial h_1\over\partial \xi_2}\right)  \eqno(184) $$
$$
\dot{\epsilon}_{22} 
= {1\over h_2}\left({\partial v\over\partial\xi_2}
   +{u\over h_1}{\partial h_2\over\partial \xi_1}\right)  \eqno(185) $$
$$
\dot{\epsilon}_{12} 
= {1 \over 2} [{h_1\over h_2}{\partial\over\partial\xi_2}\left(u\over h_1\right)
   +{h_2\over h_1}{\partial\over\partial\xi_1}\left(v\over h_2\right)]. \eqno(186) $$
Using these expressions in the dissipation $D$ (Eq. 182) and integrating by parts, 
orthogonal curvilinear forms for the stress divergence are derived:
$$
F_x = {1\over 2 }\left[{1\over h_1}{\partial\sigma_1\over \partial\xi_1} +
{1\over h_1 h_2^2}{\partial\over \partial\xi_1}\left(h_2^2\sigma_2\right) +
{2\over h_1^2 h_2}{\partial\over \partial\xi_2}\left(h_1^2\sigma_{12}\right)
\right] \eqno (187) $$
$$
F_y = {1\over 2 }\left[{1\over h_2}{\partial\sigma_1\over \partial\xi_2} -
{1\over h_1^2 h_2}{\partial\over \partial\xi_2}\left(h_1^2\sigma_2\right) +
{2\over h_1 h_2^2}{\partial\over \partial\xi_1}\left(h_2^2\sigma_{12}\right)
\right].  \eqno(188) $$

The discretization of the velocity and stress tensor components is termed the 
``bilinear discretization'' (Hunke and Dukowicz, 2002). For example, the velocity 
components are expressed in terms of the grid box vertex values and
non-dimensional coordinates ($0$ to $1$ across the grid box) as:
$$
u\left(\xi_1,\xi_2\right) 
    = u^{ne}\xi_1\xi_2 + u^{nw}\left(1-\xi_1\right)\xi_2
    + u^{sw}\left(1-\xi_1\right)\left(1-\xi_2\right)
    + u^{se}\xi_1\left(1-\xi_2\right)  \eqno(189)  $$
$$
v\left(\xi_1,\xi_2\right) 
    = v^{ne}\xi_1\xi_2 + v^{nw}\left(1-\xi_1\right)\xi_2
    + v^{sw}\left(1-\xi_1\right)\left(1-\xi_2\right)
    + v^{se}\xi_1\left(1-\xi_2\right). \eqno(190)  $$
where the four grid box velocities are referred to as ne, nw, sw, se 
for northeast, northwest, southwest and southeast respectively. Note 
that velocity is continuous across cell edges (for example, $u_{ij}^{ne} = 
u_{i+1j}^{nw}$, where $ij$ now represent grid indices). The stress tensor 
components, associated with velocity gradients through strain rates, are
discontinuous, with each cell having four corner values for stress. This 
method of discretization suppresses B-grid checkerboard solutions, 
because it is not technically B-grid, since we do not have one grid box 
center value for the stress tensor components.

The scale factors are evaluated at grid center by averaging the two
appropriate sides ($h_1=\bar{h}_1$, $h_2=\bar{h}_2$), while scale factor
spatial derivatives are simple differences in the grid side lengths 
($\partial h_1/\partial\xi_2=\Delta_2 h_1$, $\partial h_2/\partial\xi_1=\Delta_1
h_2$). Using the bilinear discretization, the strain rate terms are
evaluated as follows:

{\it Divergence}
$$
D_D^{ne}={1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{ne}-u^{nw}\right) + \Delta_1h_2u^{ne} 
+\bar{h}_1\left(v^{ne}-v^{se}\right) + \Delta_2h_1v^{ne}\right] $$  
$$
D_D^{nw}={1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{ne}-u^{nw}\right) + \Delta_1h_2u^{nw} 
+\bar{h}_1\left(v^{nw}-v^{sw}\right) + \Delta_2h_1v^{nw}\right] $$  
$$
D_D^{se}={1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{se}-u^{sw}\right) + \Delta_1h_2u^{se} 
+\bar{h}_1\left(v^{ne}-v^{se}\right) + \Delta_2h_1v^{se}\right] $$  
$$
D_D^{sw}={1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{se}-u^{sw}\right) + \Delta_1h_2u^{sw} 
+\bar{h}_1\left(v^{nw}-v^{sw}\right) + \Delta_2h_1v^{sw}\right] $$  

{\it Tension}
$$
D_T^{ne}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{ne}-u^{nw}\right) - \Delta_1h_2u^{ne} 
-\bar{h}_1\left(v^{ne}-v^{se}\right) + \Delta_2h_1v^{ne}\right] $$  
$$
D_T^{nw}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{ne}-u^{nw}\right) - \Delta_1h_2u^{nw} 
-\bar{h}_1\left(v^{nw}-v^{sw}\right) + \Delta_2h_1v^{nw}\right] $$  
$$
D_T^{se}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{se}-u^{sw}\right) - \Delta_1h_2u^{se} 
-\bar{h}_1\left(v^{ne}-v^{se}\right) + \Delta_2h_1v^{se}\right] $$  
$$
D_T^{sw}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_2\left(u^{se}-u^{sw}\right) - \Delta_1h_2u^{sw} 
-\bar{h}_1\left(v^{nw}-v^{sw}\right) + \Delta_2h_1v^{sw}\right] $$  

{\it Shearing}
$$
D_S^{ne}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_1\left(u^{ne}-u^{se}\right) - \Delta_2h_1u^{ne}
+\bar{h}_2\left(v^{ne}-v^{nw}\right) - \Delta_1h_2v^{ne}\right] $$  
$$
D_S^{nw}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_1\left(u^{nw}-u^{sw}\right) - \Delta_2h_1u^{nw}
+\bar{h}_2\left(v^{ne}-v^{nw}\right) - \Delta_1h_2v^{nw}\right] $$  
$$
D_S^{se}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_1\left(u^{ne}-u^{se}\right) - \Delta_2h_1u^{se}
+\bar{h}_2\left(v^{se}-v^{sw}\right) - \Delta_1h_2v^{se}\right] $$  
$$
D_S^{sw}= {1\over \bar{h}_1\bar{h}_2}
\left[\bar{h}_1\left(u^{nw}-u^{sw}\right) - \Delta_2h_1u^{sw}
+\bar{h}_2\left(v^{se}-v^{sw}\right) - \Delta_1h_2v^{sw}\right] $$  

For the divergence of the stress tensor, we note that for the ne
corner of grid box $ij$, there are contributions from the four
surrounding grid boxes (i.e. northeast corner of $ij$, southeast
corner of $ij+1$, northwest corner of $i+1j$ and the southwest corner
of $i+1j+1$). The terms below show these contributions to the stress
divergence for the northeast corner of grid box $ij$. Hence, for
the contributions of each grid box, the corner designations of
$\sigma$ are relative to that box. 

{\it Contribution of the $\sigma_1$ term to $F_x$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[-{\bar{h}_2\over 4}\left( 
        {1\over 3}\left(\sigma_1^{ne}+\sigma_1^{nw}\right)  
       +{1\over
6}\left(\sigma_1^{se}+\sigma_1^{sw}\right)\right)\right.\right.
-\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_1^{ne}
           +{1\over 18}\left(\sigma_1^{nw}+\sigma_1^{se}\right)
           +{1\over 36}\sigma_1^{sw}\right)\right]_{ij}  $$
$$
+\left[{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_1^{ne}+\sigma_1^{nw}\right) 
        +{1\over
6}\left(\sigma_1^{se}+\sigma_1^{sw}\right)\right)\right.
-\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_1^{nw}
           +{1\over 18}\left(\sigma_1^{sw}+\sigma_1^{ne}\right)
           +{1\over 36}\sigma_1^{se}\right)\right]_{i+1j} $$
$$
+\left[-{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_1^{se}+\sigma_1^{sw}\right) 
        +{1\over
6}\left(\sigma_1^{ne}+\sigma_1^{nw}\right)\right)\right.
-\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_1^{se}
           +{1\over 18}\left(\sigma_1^{ne}+\sigma_1^{sw}\right)
           +{1\over 36}\sigma_1^{nw}\right)\right]_{ij+1} $$
$$
+\left[{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_1^{se}+\sigma_1^{sw}\right) 
        +{1\over
6}\left(\sigma_1^{ne}+\sigma_1^{nw}\right)\right)\right.
-\left.\left.{\Delta_1h_2\over 2}
      \left({1\over 9}\sigma_1^{sw}
           +{1\over 18}\left(\sigma_1^{nw}+\sigma_1^{se}\right)
           +{1\over 36}\sigma_1^{ne}\right)\right]_{i+1j+1}\right\} $$

{\it Contribution of the $\sigma_1$ term to $F_y$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[-{\bar{h}_1\over 4}
  \left({1\over 3}\left(\sigma_1^{ne}+\sigma_1^{se}\right) 
       +{1\over
6}\left(\sigma_1^{nw}+\sigma_1^{sw}\right)\right)\right.\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_1^{ne}
           +{1\over 18}\left(\sigma_1^{nw}+\sigma_1^{se}\right)
           +{1\over 36}\sigma_1^{sw}\right)\right]_{ij} $$
$$
+\left[-{\bar{h}_1\over 4}
 \left({1\over 3}\left(\sigma_1^{nw}+\sigma_1^{sw}\right) 
      +{1\over
6}\left(\sigma_1^{ne}+\sigma_1^{se}\right)\right)\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_1^{nw}
           +{1\over 18}\left(\sigma_1^{sw}+\sigma_1^{ne}\right)
           +{1\over 36}\sigma_1^{se}\right)\right]_{i+1j} $$
$$
+\left[{\bar{h}_1\over 4}
 \left( {1\over 3}\left(\sigma_1^{ne}+\sigma_1^{se}\right) 
       +{1\over
6}\left(\sigma_1^{nw}+\sigma_1^{sw}\right)\right)\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_1^{sw}
           +{1\over 18}\left(\sigma_1^{ne}+\sigma_1^{sw}\right)
           +{1\over 36}\sigma_1^{nw}\right)\right]_{ij+1} $$
$$
+\left[{\bar{h}_1\over 4}
 \left( {1\over 3}\left(\sigma_1^{nw}+\sigma_1^{sw}\right) 
       +{1\over
6}\left(\sigma_1^{ne}+\sigma_1^{se}\right)\right)\right.
-\left.\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_1^{sw}
           +{1\over 18}\left(\sigma_1^{nw}+\sigma_1^{se}\right)
           +{1\over 36}\sigma_1^{ne}\right)\right]_{i+1j+1} \right\} $$

{\it Contribution of the $\sigma_2$ term to $F_x$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[-{\bar{h}_2\over 4}\left( 
        {1\over 3}\left(\sigma_2^{ne}+\sigma_2^{nw}\right)  
       +{1\over
6}\left(\sigma_2^{se}+\sigma_2^{sw}\right)\right)\right.\right.
+\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_2^{ne}
           +{1\over 18}\left(\sigma_2^{nw}+\sigma_2^{se}\right)
           +{1\over 36}\sigma_2^{sw}\right)\right]_{ij} $$
$$
+\left[{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_2^{ne}+\sigma_2^{nw}\right) 
        +{1\over
6}\left(\sigma_2^{se}+\sigma_2^{sw}\right)\right)\right.
+\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_2^{nw}
           +{1\over 18}\left(\sigma_2^{sw}+\sigma_2^{ne}\right)
           +{1\over 36}\sigma_2^{se}\right)\right]_{i+1j} $$
$$
+\left[-{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_2^{se}+\sigma_2^{sw}\right) 
        +{1\over
6}\left(\sigma_2^{ne}+\sigma_2^{nw}\right)\right)\right.
+\left.{\Delta_1h_2\over 2}\left({1\over 9}\sigma_2^{se}
           +{1\over 18}\left(\sigma_2^{ne}+\sigma_2^{sw}\right)
           +{1\over 36}\sigma_2^{nw}\right)\right]_{ij+1} $$
$$
+\left[{\bar{h}_2\over 4}\left( 
         {1\over 3}\left(\sigma_2^{se}+\sigma_2^{sw}\right) 
        +{1\over
6}\left(\sigma_2^{ne}+\sigma_2^{nw}\right)\right)\right.
+\left.\left.{\Delta_1h_2\over 2}
      \left({1\over 9}\sigma_2^{sw}
           +{1\over 18}\left(\sigma_2^{nw}+\sigma_2^{se}\right)
           +{1\over 36}\sigma_2^{ne}\right)\right]_{i+1j+1}\right\} $$

{\it Contribution of the $\sigma_2$ term to $F_y$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[{\bar{h}_1\over 4}
 \left( {1\over 3}\left(\sigma_2^{ne}+\sigma_2^{se}\right) 
       +{1\over
6}\left(\sigma_2^{nw}+\sigma_2^{sw}\right)\right)\right.\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_2^{ne}
           +{1\over 18}\left(\sigma_2^{nw}+\sigma_2^{se}\right)
           +{1\over 36}\sigma_2^{sw}\right)\right]_{ij} $$
$$
+\left[{\bar{h}_1\over 4}
 \left( {1\over 3}\left(\sigma_2^{nw}+\sigma_2^{sw}\right) 
       +{1\over
6}\left(\sigma_2^{ne}+\sigma_2^{se}\right)\right)\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_2^{nw}
           +{1\over 18}\left(\sigma_2^{sw}+\sigma_2^{ne}\right)
           +{1\over 36}\sigma_2^{se}\right)\right]_{i+1j} $$
$$
+\left[-{\bar{h}_1\over 4}
 \left({1\over 3}\left(\sigma_2^{ne}+\sigma_2^{se}\right) 
      +{1\over
6}\left(\sigma_2^{nw}+\sigma_2^{sw}\right)\right)\right.
-\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_2^{sw}
           +{1\over 18}\left(\sigma_2^{ne}+\sigma_2^{sw}\right)
           +{1\over 36}\sigma_2^{nw}\right)\right]_{ij+1} $$
$$
+\left[-{\bar{h}_1\over 4}
 \left({1\over 3}\left(\sigma_2^{nw}+\sigma_2^{sw}\right) 
      +{1\over
6}\left(\sigma_2^{ne}+\sigma_2^{se}\right)\right)\right.
-\left.\left.{\Delta_2h_1\over 2}
      \left({1\over 9}\sigma_2^{sw}
           +{1\over 18}\left(\sigma_2^{nw}+\sigma_2^{se}\right)
           +{1\over 36}\sigma_2^{ne}\right)\right]_{i+1j+1} \right\} $$

{\it Contribution of the $\sigma_{12}$ term to $F_x$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[-{\bar{h}_1\over 2}
  \left({1\over 3}\left(\sigma_{12}^{ne}+\sigma_{12}^{se}\right) 
 +{1\over
6}\left(\sigma_{12}^{nw}+\sigma_{12}^{sw}\right)\right)\right.\right.
+\left.{\Delta_2h_1}
      \left({1\over 9}\sigma_{12}^{ne}
           +{1\over 18}\left(\sigma_{12}^{nw}+\sigma_{12}^{se}\right)
           +{1\over 36}\sigma_{12}^{sw}\right)\right]_{ij} $$
$$
+\left[-{\bar{h}_1\over 2}
 \left({1\over 3}\left(\sigma_{12}^{nw}+\sigma_{12}^{sw}\right) 
      +{1\over
6}\left(\sigma_{12}^{ne}+\sigma_{12}^{se}\right)\right)\right.
+\left.{\Delta_2h_1}
      \left({1\over 9}\sigma_{12}^{nw}
           +{1\over 18}\left(\sigma_{12}^{sw}+\sigma_{12}^{ne}\right)
           +{1\over 36}\sigma_{12}^{se}\right)\right]_{i+1j} $$
$$
+\left[{\bar{h}_1\over 2}
 \left({1\over 3}\left(\sigma_{12}^{ne}+\sigma_{12}^{se}\right) 
      +{1\over
6}\left(\sigma_{12}^{nw}+\sigma_{12}^{sw}\right)\right)\right.
+\left.{\Delta_2h_1}
      \left({1\over 9}\sigma_{12}^{sw}
           +{1\over 18}\left(\sigma_{12}^{ne}+\sigma_{12}^{sw}\right)
           +{1\over 36}\sigma_{12}^{nw}\right)\right]_{ij+1} $$
$$
+\left[{\bar{h}_1\over 2}
 \left({1\over 3}\left(\sigma_{12}^{nw}+\sigma_{12}^{sw}\right) 
      +{1\over
6}\left(\sigma_{12}^{ne}+\sigma_{12}^{se}\right)\right)\right.
+\left.\left.{\Delta_2h_1}
      \left({1\over 9}\sigma_{12}^{sw}
           +{1\over 18}\left(\sigma_{12}^{nw}+\sigma_{12}^{se}\right)
           +{1\over 36}\sigma_{12}^{ne}\right)\right]_{i+1j+1}
\right\} $$

{\it Contribution of the $\sigma_{12}$ term to $F_y$}
$$
{1\over \bar{h}_1\bar{h}_2}\left\{ 
\left[-{\bar{h}_2\over 2}\left( 
         {1\over 3}\left(\sigma_{12}^{ne}+\sigma_{12}^{nw}\right)  
 +{1\over
6}\left(\sigma_{12}^{se}+\sigma_{12}^{sw}\right)\right)\right.\right.
+\left.{\Delta_1h_2}\left({1\over 9}\sigma_{12}^{ne}
           +{1\over 18}\left(\sigma_{12}^{nw}+\sigma_{12}^{se}\right)
           +{1\over 36}\sigma_{12}^{sw}\right)\right]_{ij} $$
$$
+\left[{\bar{h}_2\over 2}\left( 
        {1\over 3}\left(\sigma_{12}^{ne}+\sigma_{12}^{nw}\right) 
       +{1\over
6}\left(\sigma_{12}^{se}+\sigma_{12}^{sw}\right)\right)\right.
+\left.{\Delta_1h_2}\left({1\over 9}\sigma_{12}^{nw}
           +{1\over 18}\left(\sigma_{12}^{sw}+\sigma_{12}^{ne}\right)
           +{1\over 36}\sigma_{12}^{se}\right)\right]_{i+1j} $$
$$
+\left[-{\bar{h}_2\over 2}\left( 
         {1\over 3}\left(\sigma_{12}^{se}+\sigma_{12}^{sw}\right) 
        +{1\over
6}\left(\sigma_{12}^{ne}+\sigma_{12}^{nw}\right)\right)\right.
+\left.{\Delta_1h_2}\left({1\over 9}\sigma_{12}^{se}
           +{1\over 18}\left(\sigma_{12}^{ne}+\sigma_{12}^{sw}\right)
           +{1\over 36}\sigma_{12}^{nw}\right)\right]_{ij+1} $$
$$
+\left[{\bar{h}_2\over 2}\left( 
        {1\over 3}\left(\sigma_{12}^{se}+\sigma_{12}^{sw}\right) 
       +{1\over
6}\left(\sigma_{12}^{ne}+\sigma_{12}^{nw}\right)\right)\right.
+\left.\left.{\Delta_1h_2}
      \left({1\over 9}\sigma_{12}^{sw}
           +{1\over 18}\left(\sigma_{12}^{nw}+\sigma_{12}^{se}\right)
           +{1\over 36}\sigma_{12}^{ne}\right)\right]_{i+1j+1}\right\} $$

Eqns. 172,173, 178,179 and 180 are solved simultaneously over elastic time step 
$\Delta t_e = \Delta t / N_e < T_{ew} < \Delta t$ where $N_e$ is the number of 
elastic subcycle time steps. The velocities are solved for on the U-grid (i.e.
on the corners of the T-grid cells; see section 3.4). The left hand side terms 
are treated implicitly, the right hand side terms explicitly, and the rate of strain
$\dot\epsilon_{ij}$ and $\Delta$ are updated each elastic time step. The definition of
the elastic parameter $E$ in terms of the bulk viscosity, as well as the updating of 
the rate of strain tensor each elastic time step, eliminates any linearization error
associated with viscosities which are lagged over the time step.

\vskip 8pt
\noindent {\bf 4.11 Advection}

Horizontal advection in Eqs. 3-6, and 8 is evaluated after the incremental 
remapping algorithm of Lipscomb and Hunke (2004), following closely the
work of Dukowicz and Baumgardner (2000), designated DB in this section. The 
following description is adapted from the more detailed CICE documentation 
(Hunke and Lipscomb, 2004).

The remapping algorithm is second-order accurate in space, except where 
field gradients are limited to preserve monotonicity (that is, to avoid 
creating unphysical extrema or ripples in the transported fields).  One 
benefit of remapping is that it preserves the monotonicity of tracers: 
quantities that obey the advection equation $ dh/dt = \partial h / \partial t 
+ \bf{u} \cdot \bigtriangledown h = 0$, where h is a tracer.  Using the 
conservation equations 3-6 and 8, along with the definitions of ice/snow 
thickness in Eq. 21 and enthalpy in Eq. 14, one can show that thickness, 
enthalpy, and surface temperature are tracers.

The remapping algorithm can be summarized as follows: 
Given values of the ice area and tracer fields in each grid cell, construct 
linear approximations and limit the field gradients to preserve 
monotonicity. Given ice velocities at grid cell corners, identify departure
regions for the transports across each cell edge, and then divide these departure 
regions into triangles and compute the coordinates of the triangle vertices. 
Integrate these fields over the departure triangles to obtain the area, 
volume, and energy transports across each cell edge. Transfer these quantities
across cell edges and update the state variables.

Since all scalar fields are transported by the same velocity field, 
identifying departure regions is done only once per time step. The other 
three steps are repeated for each field in each thickness category. These 
steps are described below.

\vskip 10pt
\noindent {\bf 4.11.1 Reconstructing area and tracer fields}

First, using the known values of the state variables, the ice area and tracer 
fields are reconstructed in each grid cell as linear functions of $x$ and $y$. 
For each field we compute the value at the cell center (i.e., at the origin 
of a 2D Cartesian coordinate system defined for that grid cell), along with 
gradients in the $x$ and $y$ directions. The gradients are limited to preserve
monotonicity. When integrated over a grid cell, the reconstructed fields must 
have mean values equal to the known state variables denoted for example by 
$\bar{a}$ for fractional area. The mean values are not, in general, equal to 
the values at the cell center. For example, the mean ice area must equal the 
value at the centroid, which may not lie at the cell center. 

Consider first the fractional ice area. For each thickness category 
we construct a field $a({\bf r})$ whose mean is $\bar{a}$, where 
${\bf r} = (x,y)$ is the position vector relative to the 
cell center. That is, we require 
$$
\int_A a \, dA = {\bar a} \, A,  \eqno(191)    
$$
where $A=\int_A dA$ is the grid cell area.
This equation is satisfied if $a({\bf r})$ has the form
$$
a({\bf r}) = \bar{a} + \alpha \left<\nabla a\right> \cdot ({\bf r}-{\bf \bar{r}}),
  \eqno(192)    
$$
where $\left<\nabla a\right>$ is a centered estimate of the area
gradient within the cell, $\alpha$ is a limiting coefficient that
enforces monotonicity, and ${\bf \bar{r}}$ is the cell centroid:
$${\bf \bar{r}} = {1\over A} \int_A {\bf r} \, dA.$$
It follows from Eq. 192 that the ice area at the cell center 
($\bf{r} = 0$) is
$$ a_c = \bar{a} - a_x \overline{x} - a_y \overline{y}, $$
where $a_x = \alpha (\partial a / \partial x)$ and
$a_y = \alpha (\partial a / \partial y)$ are
the limited gradients in the $x$ and $y$ directions,
respectively, and the components of
${\bf \bar{r}}$, $\overline{x} = \int_A x \, dA / A$ and
$\overline{y} = \int_A y \, dA / A$. Ice thickness, snow 
thickness and enthalpy fields are treated analogously.

We preserve monotonicity by van Leer limiting. Let
$\bar{\phi}(i,j)$ denote the mean value of some field in grid cell
$(i,j)$. We first compute centered gradients of $\bar{\phi}$ in the 
$x$ and $y$ directions, then check whether these gradients 
give values of $\phi$ within cell $(i,j)$ that lie outside the
range of $\bar{\phi}$ in the cell and its eight neighbors. Let
$\bar{\phi}_{\max}$ and $\bar{\phi}_{\min}$ be the maximum and
minimum values of $\bar{\phi}$ over the cell and its neighbors,
and let $\phi_{\max}$ and $\phi_{\min}$ be the maximum and minimum
values of the reconstructed $\phi$ within the cell. Since the
reconstruction is linear, $\phi_{\max}$ and $\phi_{\min}$ are
located at cell corners. If $\phi_{\max} > \bar{\phi}_{\max}$ or
$ \phi_{\min} < \bar{\phi}_{\min}$, we multiply the unlimited
gradient by $\alpha = \min(\alpha_{\max}, \alpha_{\min})$, where
$$ \alpha_{\max} =
 (\bar{\phi}_{\max} - \bar{\phi}) / (\phi_{\max} -\bar{\phi}), $$
$$ \alpha_{\min} =
  (\bar{\phi}_{\min} - \bar{\phi}) / (\phi_{\min} -\bar{\phi}). $$
Otherwise the gradient need not be limited.

\vskip 10pt
\noindent {\bf 4.11.2 Locating departure triangles}

The locating of departure triangles is discussed in detail by DB. 
The basic idea is that the contents of a quadrilateral departure 
region are transported to the target grid cell. Consider a departure
quadrilateral north and west of the target grid cell whose motion
is towards the southeast. The neighboring grid
cells are labeled by compass directions: $NW$, $N$, $NE$, $W$, and
$E$. The four vectors point along the velocity field at the cell
corners, and the departure region is formed by joining the 
starting points of these vectors. Instead of integrating over the
entire departure region, it is convenient to compute transports across
cell edges. We identify departure regions for the north and east
edges of each cell, which are also the south and west edges of
neighboring cells.  There are 20 triangles that can contribute transports 
across the north edge of a grid cell.

This scheme was designed for rectangular grids. Grid cells in CCSM3
actually lie on the surface of a sphere and must be projected onto
a plane. Many such projections are possible. The projection used
in CCSM3 approximates
spherical grid cells as quadrilaterals in the plane tangent to the
sphere at a point inside the cell. The quadrilateral vertices are
$(N/2, E/2)$, $(-N/2,W/2)$, $(-S/2,-W/2)$, and $(S/2,-E/2)$, where
$N$, $S$, $E$, and $W$ are the lengths of the cell edges on the
spherical grid. The quadrilateral area, $(N+S)(E+W)/4$, is a good
approximation to the true spherical area. However, cell edges in
this projection are not orthogonal (i.e., they do not meet at
right angles) as on the spherical grid. This means that when
vectors are translated from cell corners to cell centers, we must
take care that the departure points in the cell-center coordinate
system lie inside the grid cell contributing the transport. Otherwise,
monotonicity may be violated, because van~Leer limiting does not
apply outside the grid cell.

Most grids cells are nearly rectangular. On the 1$^\circ$~displaced-pole 
grid used in CCSM3, the maximum angle is about 
~1$^\circ$. Vector transformations may therefore be omitted on most grids 
with little loss of accuracy. We have retained them, however, because they 
ensure exact monotonicity at little added cost.

We made one other change in the DB scheme for locating triangles. In their 
paper, departure points are defined by projecting cell corner velocities 
directly backward. That is, 
$$
\bf{x_D^\prime} = -\bf{u^\prime} \, \Delta t,  \eqno(193)    
$$
where $\bf{x}_D^\prime$ is the location of the departure point 
relative to the cell corner and the primes denote vectors defined in the 
cell-corner basis. This approximation is only first-order accurate. The 
accuracy is increased to second-order by correcting the velocity with 
a midpoint approximation before finding the departure point.

\vskip 10pt
\noindent {\bf 4.11.3 Integrating transports}

Next, we integrate the reconstructed fields over the departure
triangles to find the total transports of area, volume, and energy
across each cell edge. Area transports are easy to compute since the
area is linear in $x$ and $y$. Given a triangle with vertices
$\bf{x_i} = (x_i,y_i)$, $i\in\{1,2,3\}$, the triangle area is
$$
A_T = {1 \over 2}\left|(x_2-x_1)(y_3-y_1) -
(y_2-y_1)(x_3-x_1)\right|.  \eqno(194)    
$$
The integral $I_1$ of any linear function $f(\bf{r})$ over a
triangle is given by
$$
 I_1 = A_T f(\bf{x_0}),  \eqno(195)    
$$
where $\bf{x}_0 = (x_0,y_0)$ is the triangle midpoint,
$$
\bf{x}_0={1\over 3}\sum_{i=1}^3\bf{x}_i.  \eqno(196)    
$$
To compute the area transport, we evaluate the area at the midpoint,
$$
a(\bf{x}_0)  = a_c + a_x x_0 + a_y y_0,  \eqno(197)    
$$
and multiply by $A_T$. By convention, northward and eastward
transports are positive, while southward and westward transports are
negative.

Eq. 197 cannot be used for volume transports, because the reconstructed volumes 
are quadratic functions of position. (They are products of two linear 
functions, area and thickness.) The integral of a quadratic polynomial 
over a triangle requires function evaluations at three points,
$$
 I_2 = {A_T \over 3}\sum_{i=1}^3 f\left({\bf x}^\prime_i\right),  \eqno(198)    
$$
where $\bf{x}_i^\prime = (\bf{x}_0+\bf{x}_i)/2$ are points lying halfway 
between the midpoint and the three vertices. Eq. 198 does not work for ice 
and snow energies, which are cubic functions---products of area, thickness, 
and enthalpy. Integrals of a cubic polynomial over a triangle can be 
evaluated using a four-point formula:
$$
 I_3 = A_T \left[ -{9 \over 16} f(\bf{x}_0) +
              {25 \over 48} \sum_{i=1}^3 f(\bf{x}_i^{\prime\prime})\right]  \eqno(199)    
$$
where $\bf{x_i}^{\prime\prime}=(3 \bf{x}_0 + 2
\bf{x}_i)/5$.

\vfill
\eject
\vskip 10pt
\noindent {\bf 4.11.4. Updating state variables}

Finally, we use these transports to compute new values of the state variables
in each ice category and grid cell. The new fractional ice areas $a_{in}^\prime(i,j)$ 
are given by
$$ 
a_{in}^\prime(i,j) = a_{in}(i,j) +
              {{F_{Ea}(i-1,j) - F_{Ea}(i,j)
                  + F_{Na}(i,j-1) - F_{Na}(i,j)} \over
                   {A(i,j)}}  \eqno(200)    
$$
where $F_{Ea}(i,j)$ and $F_{Na}(i,j)$ are area transports across the
east and north edges, respectively, of cell $(i,j)$, and $A(i,j)$
is the grid cell area. All transports added to one cell are subtracted 
from a neighboring cell; thus Eq. 200 conserves total 
ice area.

The new ice volumes and internal energies are computed analogously. New 
thicknesses are given by the ratio of volume to area, and enthalpies
by the ratio of energy to volume. Tracer monotonicity is ensured because
$$ h^\prime = {\int_A a \, h \, dA \over \int_A a \, dA}, $$
$$ q^\prime  = {\int_A a \, h \, q\,dA \over \int_A a \, h \ dA}, $$
where $h^\prime$ and $q^\prime$ are the new-time thickness and
enthalpy, given by integrating the old-time ice area, volume, and
energy over a Lagrangian departure region with area $A$. That is,
the new-time thickness and enthalpy are weighted averages over
old-time values, with non-negative weights $a$ and $ah$. Thus the
new-time values must lie between the maximum and minimum of the
old-time values.

\vskip 8pt
\noindent {\bf 4.12 Mechanical Redistribution}

Mechanical redistribution of ice thickness due to rafting and ridging
processes is treated in this section. Specifically, the redistribution
function $R$ of Eq. 1 ($\psi$ in Thorndike et al. 1975)
is parameterized, and then used to evaluate the $S_M$ source/sink terms 
in Eqs. 3-6,8. For example, the integral of $R$
over the $n^{th}$ category gives the ice fraction source
$$
    S_{MAn}=\int_{h_n^*}^{h_{n+1}^*} R(h,g,\bf{u}) dh .  \eqno(201)   
$$

The redistribution function $R$ depends on the the ice thickness $h$, 
the distribution function $g(h)$, and the velocity field $\bf{u}$,
specifically the invariants of the strain rate tensor (see Section 4.10,
Eqs 159,165,166). These invariants are the divergence $\dot \epsilon_I$ 
and the shear $\dot \epsilon_{II}$. These two invariants are also used in 
terms of magnitude $\vert \dot \epsilon \vert=( \dot \epsilon_I^2+\dot \epsilon_{II}^2)^{1/2}$ 
and strain rate angle $\theta=\tan^{-1} (\dot \epsilon_{II}/\dot \epsilon_I)$.
(Note that $\theta=0^{\circ}$ refers to pure divergence,
$\theta=45^{\circ}$ uniaxial extension, $\theta=90^{\circ}$ pure shear,
$\theta=135^{\circ}$ to uniaxial compression, and $\theta=180^{\circ}$ pure convergence). 

Two strong constraints on the redistribution follow from the first two 
moments of the distribution equation, corresponding to area and volume 
conservation:
$$
\int_0^{h_{max}} R dh = \dot\epsilon_I  \eqno(202)   
$$
$$
\int_0^{h_{max}} h R dh = 0.  \eqno(203)   
$$
The first follows from the conservation of total area, since the import or
export of area ($\dot\epsilon_I$) must be balanced by a redistribution. 
The second follows from the requirement that mechanical redistribution
cannot change the ice volume.

The parameterization of lead opening and mechanical redistribution
follows from the theoretical formulation of Thorndike et al. (1975):
$$
  R=\vert \dot\epsilon \vert[\alpha_0(\theta) \delta(h)+\alpha_r(\theta) w_r(h,g)],  \eqno(204)   
$$
where $\delta(h)$(the delta function) is the opening mode, and $w_r(h,g)$ 
is the ridging mode. 
Note that the conservation of area requires:
$$
\int_{0}^{\infty}w_r(h,g)dh=-1.  \eqno(205)   
$$
The coefficients $\vert \dot\epsilon \vert \alpha_0(\theta)$ and $\vert \dot\epsilon \vert
\alpha_r(\theta)$ are known as the lead opening and closing rates, respectively, and they 
are related such that their difference equals the ice divergence, $\vert \dot\epsilon \vert 
\alpha_0(\theta)-\vert \dot\epsilon \vert \alpha_r(\theta)=\dot\epsilon_I$.  

Two aspects of the mechanical redistribution must be considered: the relation
between the distribution function and the deformation/compressive strength used 
in the ice dynamics, and the redistribution source terms in the conservation equations.

For the first, we follow an energetics argument by Rothrock (1975). The deformational
work done on the ice is equated to known sinks of energy in ridge building:
$$
  \sigma_I \dot \epsilon_{I} + \sigma_{II} \dot \epsilon_{II} = R_{pot} + R_{fric}  \eqno(206)    
$$
where $R_{pot}$ is the rate of mechanical production of gravitational potential energy
per unit area, and $R_{fric}$ is the rate of frictional energy loss per unit area. We 
write
$$
  \sigma_I \dot \epsilon_{I} + \sigma_{II} \dot \epsilon_{II} = 
(1 + {{R_{fric}} \over {R_{pot}}})R_{pot} = ZR_{pot}  \eqno(207)    
$$
where $Z$ is the ratio of total energy dissipated to potential energy gain. Ice thickness
$h$ has potential energy relative to sea level of $P_e = C_{pe} h^2$, where 
$C_{pe} = {1 \over 2} {{\rho_i} \over {\rho_o}}(\rho_o - \rho_i) g_e $ 
Integrating over the entire thickness distribution results in the total potential energy
$$
P_e = C_{pe} \int_0^\infty h^2 g(h) dh  \eqno(208)    
$$
The rate of gain of $P_e$ is
$$
{dP_e \over dt} = C_{pe} \int_0^\infty h^2 {{dg(h)} \over {dt}} dh .  \eqno(209)    
$$
From the distribution equation, using $\bigtriangledown\cdot \bf{u}g =
g\bigtriangledown\cdot \bf{u} + \bf{u}\cdot\bigtriangledown g$:
$$ 
    {{dg} \over {dt}} = - {{\partial} \over {\partial h}} ({\dot{h}}g) 
    + L(h,g) - g\bigtriangledown\cdot \bf{u} + R .  \eqno(210)    
$$
Thus the total rate of change of potential energy is
$$ 
       {{dP_e} \over {dt}} = 
- \int_0^\infty h^2 {{\partial} \over {\partial h}} ({\dot{h}}g)dh
+ \int_0^\infty h^2 L(h,g)dh - P_e\bigtriangledown\cdot \bf{u} 
+ \int_0^\infty h^2 R dh .  \eqno(211)    
$$
The first two terms on the right hand side refer to gain/loss due to
thermodynamic and thickness transport processes, the third to large
scale divergence, while the last term is due to mechanical redistribution. 
Using Eq. 204 for $R$ we have
$$ 
R_{pot} = \vert \dot\epsilon \vert \alpha_r(\theta) C_{pe} 
\int_0^\infty h^2 w_r(h,g) dh .  \eqno(212)    
$$
We can therefore define 
$$ 
P = Z C_{pe} \int_0^\infty h^2 w_r(h,g) dh  \eqno(213)    
$$
so that
$$
   \sigma_I \dot \epsilon_{I} + \sigma_{II} \dot \epsilon_{II} = \vert \dot\epsilon \vert 
   \alpha_r(\theta) P  \eqno(214)    
$$
where $P$ is the compressive strength used in the dynamics (Eq. 163).

Further, from the viscous/plastic relations of $\sigma_I$ and $\sigma_{II}$ in terms
of $\dot \epsilon_{I}$ and $\dot \epsilon_{II}$, and the bulk/shear viscosity 
definitions (see Section 4.10, Eqs. 161,162,167), we have:
$$
   \sigma_I \dot \epsilon_{I} + \sigma_{II} \dot \epsilon_{II} = 
   \zeta \dot \epsilon_{I}^2 - {P \over 2}\dot \epsilon_{I} + 
   \eta \dot \epsilon_{II}^2   \eqno(215)    
$$
or
$$
   \sigma_I \dot \epsilon_{I} + \sigma_{II} \dot \epsilon_{II} = P \vert \dot\epsilon \vert
   \{ {1 \over 2\vert \dot\epsilon \vert}(\Delta - \dot \epsilon_{I})  \eqno(216)    
$$
so that the ridging mode $\alpha_r(\theta)$ in terms of the strain rate angle is
$$
   \alpha_r(\theta) = -{1 \over 2}cos(\theta) + {1 \over 2} \sqrt{cos^2(\theta)
   + {{sin^2(\theta)} \over {e^2}}}  \eqno(217)    
$$
$$
  \vert \dot\epsilon \vert \alpha_r(\theta)={1 \over 2} (\Delta-\dot \epsilon_I).
  \eqno(218)    
$$
which is the result for the elliptical yield curve with aspect ratio $e$ found
by Hibler (1980). 

Flato and Hibler (1995) separated the expression in Eq. 218
into the sum of two terms representing the energy dissipation from ridge
building by shear and convergence:
$$
\vert \dot\epsilon \vert \alpha_r(\theta)=C_s {1 \over 2} (\Delta-\mathopen|\dot
\epsilon_I \mathclose|)-{\rm min}(\dot \epsilon_I,0),  \eqno(219)    
$$
with the factor $C_s$ added so the shearing component can be altered
by varying $C_s$ between 0 (all energy dissipation by shear is lost to
sliding) and 1 (all energy dissipation by shear is used to build ridges).  
Eqs. 218 and 219 are equivalent when $C_s=1$. The experiments of Flato 
and Hibler (1995) indicate that $C_s=0.5$
is appropriate to produce the concentrations of ridged ice observed in
the Arctic. However, Bitz et al. (2001) found that a coupled model tended
to predict too much ridged ice with $C_s=0.5$.  Bitz et al. (2001) tested 
the model with $C_s=0$ and found better agreement with observations.
Unfortunately, the parameter $C_s$ depends on $e$ and $P$, and none of
these values is well constrained by observations.  It is some
consolation that Bitz et al. (2001) found that $C_s$, and hence the precise
concentration of ridged ice, has relatively little affect on the {\it
climate} of the Arctic. Our standard model uses a compromise value
of $C_s=0.225$. 

The ridging mode is the sum of two distributions describing the ice
participating in ridging $a(h)$ and the newly ridged ice
$n(h)$, normalized to conserve area and volume:
$$
  w_r(h)= { -a(h)+ n(h) \over \int^\infty_0 [a(h)-n(h)] dh}.  \eqno(220)    
$$
The ice participating in ridging is found by weighting $g(h)$ by a
function $b(h)$ that is designed to make thinner ice more likely to
ridge than thicker ice.  The newly ridged ice is found by integrating
$a(h)$ times the redistribution function $\gamma(h',h)$ over the range
of ice thicknesses $h'$ that can contribute to form newly ridged ice
of thickness $h$. Hence,
$$
  a(h) =  b(h)g(h) 
  n(h) =  \int^h_0 \gamma(h',h) b(h')g(h') dh'.   \eqno(221)    
$$
Thorndike et al. (1975) argued that a plausible $b(h)$ might depend
linearly on the cumulative thickness distribution $G(h)=\int_0^h g(h')
dh'$ according to
$$
b(G)= {2 \over G^*} \left[1-{G(h) \over G^*}\right]  \eqno(222)    
$$
when $G \le G^*$; otherwise $0$ for $G > G^*$, for which
$G^*$ is the limiting fraction below which all ridging occurs
and is assumed to be 15\% as in Thorndike et al. (1975).  The
redistribution process is parameterized according to Hibler (1980), 
who constructed a rule for deriving $n(h)$ from $a(h)$ based on 
observations that constrain ice of thickness $h$ participating in ridging 
to be linearly distributed between thicknesses $2h$ and $2\sqrt{Kh}$:
$$
\gamma(h',h)= {1 \over 2(K-h')}  \eqno(223)    
$$
when $2h' \le h \le 2\sqrt{Kh'}$; otherwise $0$. 

The mechanical redistribution function $R_n$ is the integral of the continuous 
function (see Eq. 204) over the thickness limits of each category: 
$$
R_n=\int_{h_n^\ast}^{h_{n+1}^\ast} R dh = \delta(h) 
  \left[\dot \epsilon_I +  \vert \dot\epsilon \vert \alpha_r(\theta) \right]
  + \vert \dot\epsilon \vert \alpha_r(\theta) W_n,   \eqno(224)    
$$
where
$$
W_n=\int_{h_n^\ast}^{h_{n+1}^\ast} w_r(h) dh.   \eqno(225)    
$$
The $W_n$ factors can be separated into two components, participation
$W_{an}$ (loss) and redistribution $W_{nn}$ (gain),
$$
W_n= -W_{an} + W_{nn},   \eqno(226)    
$$
such that % changes begin here
$$
\eqalign{
W_{an}&={1 \over \omega} \int_{h_n^\ast}^{h_{n+1}^\ast} b(h)g(h)dh \cr
W_{nn}&={1 \over \omega} \int_{h_n^\ast}^{h_{n+1}^\ast} \int_0^{h_{n+1}^\ast} 
 \gamma(h',h) b(h')g(h')dh'dh  \cr }   \eqno(227)    
$$ 
where $\omega$ normalizes $W_n$ such that $\sum_{n=0}^{N}R_n=\dot
\epsilon_I$. After substituting $g(h)dh=dG$ into the equation for
$W_{an}$, we find
$$
  W_{an} = {1 \over \omega} \int_{{\rm min}(G^*,G_n)}
  ^{{\rm min}(G^*,G_{n+1})}b(G)dG,   \eqno(228)    
$$ 
where $G_n= \sum_{p=0}^n A_p$, with $G_{-1}=0$. Finally we express
$W_{an}$ as a function of the auxiliary function $Y_n$:
$$
 W_{an} = {1 \over \omega}(Y_n-Y_{n+1}), \ \ {\rm where } \ \ 
 Y_n=
       [1-{G_n \over G^*}]^2   \eqno(229)    
$$ 
when $G_n \le G^*$; otherwise $0$ for $G_n > G^*$. 
In Eq. 227, $W_{nn}$ is evaluated by first changing the order of
integration and then expanding the outer integral into a sum of
integrals over the categories:
$$
    W_{nn}= {1 \over \omega} \sum^{n+1}_{p=0} \int_{h_p^\ast}^{h_{p+1}^\ast} 
      \int_{h_n^\ast}^{h_{n+1}^\ast} \gamma(h',h) dh ~ b(h')g(h')dh'.   \eqno(230)    
$$
Taking $\int_{h_n^\ast} 
  ^{h_{n+1}^\ast} \gamma(h',h) dh $ outside of the remaining
integral we have
$$
    W_{nn}= {1 \over \omega} \sum^{n+1}_{p=0}   \Gamma_{pn+1}
     \int_{h_n^\ast}^{h_{n+1}^\ast} b(h')g(h')dh' = 
    {1 \over \omega} \sum^{n+1}_{p=0} W_{ap} \Gamma_{pn+1}.   \eqno(231)    
$$
The discrete distributor $\Gamma_{pn+1}$ is computed as in the method of
Hibler (1980) from
$$
\Gamma_{pn+1} = \int_{h_n^\ast}^{h_{n+1}^\ast} \gamma(h_p,h) dh =
{{\rm min}(2\sqrt{Kh_p},h_{n+1}^\ast)-{\rm max}(2h_p,h_n^\ast) \over
  2(K-h_p) },   \eqno(232)    
$$
when $2h_p<h_n^\ast$ or $2\sqrt{Kh_p}>h_{n+1}^\ast$; otherwise
$\Gamma_{pn+1}=0$. 
In general, ice from any category may participate in ridging provided 
that the cumulative distribution of ice up to that category is less 
than $G^*$. $S_{MAn}$ (Eq. 201) can now be expressed as:
$$
\eqalign{
  S_{MAn}  &= \vert \dot\epsilon \vert \alpha_r(\theta)
    \left[ - W_{an} + \sum_{p=1}^{n+1} 
    W_{ap}\Gamma_{pn+1} \right] . \cr }   \eqno(233)    
$$
Because the volume of the ice and snow is influenced by ridging, the
volume of the ice and snow must be redistributed as well as the ice
concentration. The ice participating in ridging reduces the volume in 
category $n$ proportional to $h_n W_{an}$. The ice that ridges into category 
$n$ increases the volume proportional to $\sum_{p=1}^{n+1} \tilde h_{pn+1}
W_{an}\Gamma_{pn+1}$, which is summed over the $p$ categories that ridge
into category $n$. The thickness of this newly ridged ice $\tilde
h_{pn+1}$ is uniquely determined by conservation of volume and the
discrete distributor. For simplicity, we assume the snow thickness
redistribution process is independent of the category that receives
the newly ridged ice and snow.  Hence, the snow thickness would be the
same on top of each ridged ice category. The rate of change of the ice 
and snow volumes due to mechanical redistribution is
$$
\eqalign{
  S_{MVn}  &= \vert \dot\epsilon \vert \alpha_r(\theta)
    \left[ - h_n W_{an} + \sum_{p=1}^{n+1} \tilde h_{pn+1}
    W_{ap}\Gamma_{pn+1} \right] \cr
  S_{MVsn} &= \vert \dot\epsilon \vert \alpha_r(\theta)
   \left[ -h_{sn} W_{an} + 
   \sum_{p=1}^{n+1} \tilde h_{spn+1} W_{ap}\Gamma_{pn+1} \right], }   \eqno(234)    
$$
where
$$
 \tilde h_{pn+1} ={1 \over 2}  \left[{\rm max}\left(2h_p,h_n^\ast\right)
     +{\rm min}\left(2\sqrt{Kh_p},h_{n+1}^\ast\right)\right] 
 \tilde h_{spn+1}=h_{sp} { h_p+\sqrt{Kh_p}\over h_p}.   \eqno(235)    
$$
Conservation of energy requires a redistribution of internal energy as
well.  We assume that mechanical redistribution does not mix energy
vertically.  Again it may be helpful to think of the vertical
dimension broken into a fixed number of layers, where energy is
redistributed layer by layer. Hence, redistribution transfers heat
only from the upper layer of one category to another and so on.
Consistent with the redistribution of ice and snow volume, we find
$$
S_{MEn}  = \vert \dot\epsilon \vert \alpha_r(\theta)
 \left[ - q_n h_n W_{an} + 
   \sum_{p=1}^{n+1} q_p
   \tilde h_{pn+1} W_{ap}\Gamma_{pn+1}\right]  \eqno(236)
$$
$$
S_{MESn} = \vert \dot\epsilon \vert \alpha_r(\theta)
 \left[ - q_s h_{sn} W_{an} + 
 \sum_{p=1}^{n+1} q_{sp} \tilde h_{spn+1} W_{ap}
   \Gamma_{pn+1}\right].   \eqno(237)    
$$
The ridged ice will have the same vertical temperature profile as the 
ice which has participated in ridging.

Finally, the surface temperature is also affected by mechanical
redistribution:
$$
 S_{MTsn} = \vert \dot\epsilon \vert \alpha_r(\theta) \left[-T_{sn}W_{an}+\sum_{p=1}^{n+1}
   T_{sp}W_{ap}\Gamma_{pn+1}\right].    \eqno(238)    
$$
With the redistribution of surface temperature, the ridged ice will have 
the same surface temperature as the ice that participates in ridging.

Using the definitions above, the mechanical pressure P (see Eq. 163) is:
$$
  P=Z C_{pe} \sum_{n=1}^N \left[ -h_n^2 W_{an}+ \sum_{p=1}^{n+1} \tilde
  h^2_{pn+1} W_{ap}\Gamma_{pn+1}\right].   \eqno(239)    
$$
Based on the work of Hopkins and Hibler (1991) and Flato and Hibler (1995),
we let $Z=17$.

\vfill
\eject
\vskip 8pt
\noindent {\bf 4.13 Output to the Coupler}

Aggregate states and fluxes over the ice distribution are computed
for exchange with the coupler and for history output. The general 
aggregate equation for an arbitrary field ($\{X_n\}, n=1,...N$) is;
$$  
               X = {1 \over A} {\sum_{n=1}^N}X_n A_n     \eqno(240)   
$$
where $A = \Sigma_{n=1}^N A_n$ (note that some history fields do
not normalize by $1/A$; see section 6).

Table 3 lists the states and fluxes which are sent by the sea ice model 
to the coupler (see section 2.3). The ice area $A$ follows from Eqs. 2 and 3.
The surface temperature $T_s$ comes from the vertical heat conduction
solutions in sections 4.6.1-5 . The albedos are from Eq. 35 (direct
and diffuse albedos are identical; see section 4.4). 

The latent and sensible heat fluxes $F_{LH}$ and $F_{SH}$ are from 
Eqs. 44 and 43 respectively. The upwelling flux $F_{LWUP}$ is from 
Eq. 40. The evaporative flux $F_{EVAP}$ is $F_{LH}$ without the
latent heat of sublimation $L_s$ in Eq. 44. 

The atmosphere/ice stresses $\tau_{ax}$ and $\tau_{ay}$ are given by
Eqs. 41 and 42 respectively. These stresses are computed on the T-grid 
in the ice model, and are rotated from the displaced pole grid onto the 
geographical latitude/longitude directions before being sent to the 
coupler. Similarly, the ocean/ice stresses $\tau_{ox}$ and $\tau_{oy}$,
which are computed on the U-grid from Eq. 157, are first bilinearly
interpolated to the T-grid, and then reprojected onto the geographical 
latitude/longitude directions. In this manner, all vector fields 
exchanged with the coupler are on the T-grid and projected onto 
geographical latitude/longitude directions.

Specifically, the following U-grid to T-grid interpolation is required. Note 
that ($u,v$) represents a vector field for which the subscripts ``$g$'' and ``$dp$''
refers to geographic and displaced pole grids respectively, and the superscripts 
``$t$''and ``$u$'' to the T-grid and U-grid respectively. 
$$
\eqalign{
      (u_{dp}^t)_{ij}  = {1 \over 4}(A^u_{ij}(u_{dp}^u)_{ij} + A^u_{i-1j}(u_{dp}^u)_{i-1j}
      + A^u_{ij-1}(u_{dp}^u)_{ij-1} + A^u_{i-1j-1}(u_{dp}^u)_{i-1j-1}) / A^t_{ij} \cr
      (v_{dp}^t)_{ij}  = {1 \over 4}(A^u_{ij}(v_{dp}^u)_{ij} + A^u_{i-1j}(v_{dp}^u)_{i-1j}
      + A^u_{ij-1}(v_{dp}^u)_{ij-1} + A^u_{i-1j-1}(v_{dp}^u)_{i-1j-1}) / A^t_{ij} \cr
}     \eqno(241)   
$$
where $A^u_{ij}$ is the U-grid box area, and $A^t_{ij}$ is the T-grid 
box area.

Then, the vector components ($u_{dp}^t,v_{dp}^t$) are rotated back to the geographic
grid orientation by:
$$
\eqalign{
      (u_g^t)_{ij}  &= (u_{dp}^t)_{ij} cos(\chi^t_{ij}) - (v_{dp}^t)_{ij} sin(\chi^t_{ij}) \cr
      (v_g^t)_{ij}  &= (u_{dp}^t)_{ij} sin(\chi^t_{ij}) + (v_{dp}^t)_{ij} cos(\chi^t_{ij}) \cr
}     \eqno(242)   
$$
where ($ij$) are the longitude/latitude indices of the appropriate grid. 

The shortwave flux transmitted to the ocean is:
$$
F_{SWon}  = I_{0vs} e^{-\kappa_{vs} h_n} +I_{0ni} e^{-\kappa_{ni} h_n}    \eqno(243)   
$$
(see Eq. 75 and following).

The heat flux to the ocean $F_{Qio}$ (Eq. 136) represents that heat flux available 
$F_{Qoi} < 0$ that is actually used in ice melt. The heat used by the ice model
includes basal and lateral melt.

For the water flux $F_{Wo}$, recall that
three types of ice formation are distinguished: {\bf{frazil}} (which forms 
directly in the ocean surface layer), {\bf{congelation}} (which forms at ice base), 
and {\bf{snow-ice}} (which forms by flooding of snow covered ice). Frazil ice is 
formed in the ocean model as previously described (see Section 4.3). The amount of
frazil ice given to the ice model is implied in the freezing/melting potential
($F_{Qoi}>0$). Therefore, there is no explicit water flux from the ocean to the 
ice model in this case. However, when the ice melts, the melt water is 
passed back to the ocean. Congelation ice forms in the ice model at ice base. 
Therefore, the water exchange with the ocean must include both formation and melt. 
Snow-ice is an internal transformation in the ice model itself, and need not be 
included in the water flux to the ocean.

The water flux to the ocean for category $n$ is
$$  
F_{Won} = A_n F_{RN} + \{ A_n(-\rho_i \delta h_t -\rho_i \delta h_b
-\rho_s \delta h_s) + (R_{side}(\rho_i V+ \rho_s V_s))_n \} / \Delta t    \eqno(244)   
$$  
where $A_nF_{RN}$ is the flux due to rain on ice assumed to drain directly to the 
ocean, $\delta h_t$ is the change in the top ice thickness due to surface melt and 
sublimation/condensation, $\delta h_b$ is the change in the basal ice thickness 
due to congelation ice formation and melt, $\delta h_s$ is the change in surface 
snow depth due to melting and evaporation, and $R_{side}(\rho_i V+ \rho_s V_s)$ 
is the ice and snow amount melted by lateral thermodynamic processes (see 
section 4.8).

Presently, frazil ice formation occurs within the ocean, and therefore the amount
of salinity in frazil ice must be known to the ocean. This value is the reference
salinity $S_i$, which then becomes the constant salinity that all ice must have for
salinity conservation in ocean-ice exchange. This requires that non-frazil ice
formation and all ice melt have salinity exchanges with the ocean (see Eqs.
123-126). The resulting ocean-ice salinity flux for category $n$ is
$$  
F_{S_{on}} = \{ \rho_i S_i \delta V_{ice~melt} -\rho_i S_i \delta V_{cong~ice} 
-\rho_i S_i \delta V_{sub~cond} -\rho_i S_i \delta V_{snow-ice} \}_n / \Delta t \eqno(245)   
$$  
where the melt of total ice volume $\delta V_{ice~melt}$, congelation ice formation 
$\delta V_{cong~ice}$, and snow-ice formation $\delta V_{snow-ice}$ are each positive, 
and the atmosphere-ice sublimation $\delta V_{sub~cond}$ is negative for sublimation 
and positive for condensation. It is admitted that salinity exchange with the ocean
for sublimation/condensation and for snow-ice formation does not seem very physical,
but it is required by the present version.

Diagnostic fields of 2 m reference temperature $T_{ref}$ and specific humidity 
$Q_{ref}$ follow from Eqs. 67-68 and 70 respectively. The total absorbed shortwave
flux in the ice/ocean system derives from Eq. 39.

\vfill
\eject
\noindent {\bf 5. Active Ice Only (AIO) Framework}

The ice model can be run through the coupler but with other components prescribed in
a framework called Active Ice Only (AIO). In AIO the ice model communicates via 
the coupler with other components in the same manner as it would in a fully coupled 
run. The user need only ensure that the other components are data models (i.e. prescribed), 
and that they provide the sea ice model with reasonable values of all fields received from 
the coupler (see Table 2). Because of the way ice and ocean are coupled, the ice-ocean 
heat exchange in the data ocean model is independent of ice state. Therefore
the sea ice model has an option to run with a simple ocean mixed layer model that
is part of the sea ice component.

The mixed layer prognostic variable is ocean temperature $T_o$, determined by the 
thermodynamic equation:
$$
   \rho_o c_o h_o {{\partial T_o} \over {\partial t}} = F_{SW} + F_{LW} 
        + F_{SH} + F_{LH} + F_{Qio} - F_{Qo}.   \eqno(246)
$$
where $\rho_o, c_o, h_o$ are the ocean mass density, heat capacity and mixed layer 
depth respectively, $F_{SW}$ is the absorbed shortwave flux, $F_{LW}$ is the 
absorbed longwave flux: $F_{LW} = F_{LWDN} - F_{LWUP}$, $F_{SH}$,$F_{LH}$ are
the sensible and latent heat fluxes between ocean and atmosphere respectively,
$F_{Qio}$ is the available mixed layer heat flux used by the sea ice model as 
shown in section 4.8, and $F_{Qo}$ is the ocean heat flux which represents 
seasonal mixed layer heat storage/release and lateral oceanic heat transport. Note that 
the fluxes include contributions from both open water and sea ice. The required input 
(presently monthly mean) data are the mixed layer depth, salinity (presently not used), 
ocean currents, sea surface slopes and ocean heat flux, as shown in Table 6. The 
ocean currents and tilt are used in the sea ice dynamic calculation. 

All of these data are supplied on a monthly basis from either an observational 
and/or a model source. For model data, the mean monthly surface temperatures 
and surface heat fluxes from a model run can be used. $F_{Qo}$ is computed from 
the monthly data as:
$$
   F_{Qo}^k = -\rho_o c_o h_o^k \left( {{T_o^{k+1}-T_o^{k-1}} \over {2 \Delta t_k}} 
               \right) +
               F_{SW}^k + F_{LW}^k + F_{SH}^k + F_{LH}^k + F_{Qio}^k \eqno(247)
$$
for the $k^{th}$ month ($k$ = 1, 2, ..., 12), and the required data are shown, with 
$\Delta t_k = $ mean month time. The monthly data (if necessary) must be spatially 
interpolated to the sea ice model grid, and be formatted as an external netCDF file,
containing the monthly fields shown in Table 6. Once read in, these monthly data 
are linearly interpolated in time to the sea ice model time step (therefore allowing 
for seasonal variation only).

\vfill
\eject
\vbox{
\hang \noindent {{\it Table 6.} 
Ocean Fields Required for Mixed Layer}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Symbol & Description & Units \cr
\noalign{\hrule height 1pt}
      $h_o$      &  mixed layer depth            &    m               \cr
      $S_o$      &  salinity                     &    ppt             \cr
      $u_o$      &  x direction surface velocity &    m s$^{-1}$      \cr
      $v_o$      &  y direction surface velocity &    m s$^{-1}$      \cr
$(\nabla H_o)_x$ &  x direction surface slope    &    m m$^{-1}$      \cr
$(\nabla H_o)_y$ &  y direciton surface slope    &    m m$^{-1}$      \cr
      $F_{Qo}$   &  ocean heat flux              &    W m$^{-2}$      \cr
\noalign{\hrule height 1pt}
}}}

For a time step $m$ with initial temperature $ T_o^m$, the 
mixed layer temperature forecast equation is:
$$
    T_o^{'m+1} = T_o^m + \Delta t \left( F_{SW}^m + F_{LW}^m + F_{SH}^m 
               + F_{LH}^m + F_{Qio}^m \right) / \rho_o c_o h_o   \eqno(248)
$$
where we first evaluate the exchange between the mixed layer and the atmosphere/ice 
above, in order to limit possible loss of mixed layer heat to the deep heat source when
the mixed layer temperature is at freezing (see below). $F_{SW}^m$ is computed from:
$$
      F_{SW}^m = F_{SWDN}^m (1-\alpha_o)(1-A^m) + F_{SW_o}^m A^m  \eqno(249)
$$
where $\alpha_o$ is a constant ocean surface albedo, $F_{SW_o}^m$ is the shortwave
flux that penetrates the ice to be absorbed in the underlying ocean, and $A^m$ is
the sea ice fraction. The fluxes $F_{LWUP}^m$, $F_{SH}^m$ and $F_{LH}^m$ are computed 
over the open ocean using the ocean mixed layer temperature and surface properties, 
and then weighted by the open ocean fraction. 

If $T_o^{'m+1} < T_{of}$, where $T_{of}$ is the freezing temperature of the ocean, 
and $F_{Qo} > 0$ (implying heat loss from the mixed layer), then the heat exchange
is limited such that $F_{Qo}=0$. In other cases, deep exchange is evaluated by:
$$
      T_o^{m+1} = T_o^{'m+1} - F_{Qo}^{*} \Delta t /\rho_o c_o h_o . \eqno(250)
$$
Frazil ice heat flux ($>$0) or melt potential ($<$0) is then evaluated as:
$$
      F_{Qoi} = \rho_o c_o h_o {{(T_{of} - T_o^{m+1})} \over {\Delta t}}  \eqno(251)
$$
as in the ocean model (see Section 4.3). If $T_o^{m+1} < T_{of}, 
T_o^{m+1} = T_{of}$, to ensure the mixed layer temperature always remains 
above freezing.

\vfill
\eject
\noindent {\bf 6. Output to History Files}

We present the history fields typically written to history file(s) during
an ice model integration in order to relate them to named variables used
in this technical note. Averaged quantities are written to the history file 
at the end of the time stepping loop.

There are two groups of fields written to the history file: time-invariant and
time dependent. The former are fields related to the horizontal grid, while the 
latter are prognostic and diagnostic fields from an ice model integration.
Table 7 presents a list of all history fields; the first fields from tmask
to ANGLET are time-invariant, while the rest are prognostic or diagnostic.

The field names and the equivalent variable name used in 
this document (in parentheses) are shown. Fields received from 
or sent to the coupler have (cpl) in their description. Note that some 
of the fields have non-standard units and names not necessarily consistent 
with this document. 

In some cases, two values of a flux are written to the history file: the
value sent to the coupler and a value not divided by the
ice area $A$ (see Eq. 240). The value sent to the coupler has been divided by ice
area, since the coupler multiplies by ice area, while the value used by
the ice model does not.

History file format is netCDF. The history file frequency
(i.e. how often written) can be either instantaneous (i.e. every time step),
daily, monthly, or annual. A sequence of history files is written during
model execution at the desired frequency, and all but the instantaneous files 
are time averaged over the interval between writes (see the CSIM User's Guide 
on the CCSM web page for more information). There is one exception to time 
averaging: the normalized principal stress components (sig1 and sig2) are 
always instantaneous because time average stress states at geographically 
fixed points are not physically meaningful.

There are four tendency fields included in the history file: two ice volume
and two ice area tendencies [$(\partial V / \partial t)_T$, 
$(\partial V / \partial t)_D$] and [$(\partial A / \partial t)_T$), 
$(\partial A / \partial t)_D$] respectively. These tendencies are purely
diagnostic, and distinguish the effect of thermodynamic processes (vertical
and lateral, designated by subscript T) from dynamic processes (advection 
and rafting/ridging, designated by subscript D). For example, the ice volume 
change due to thermodynamic processes is:
$$
(\partial V / \partial t)_T = [ V(after~vertical,~lateral~thermodynamics) - V(before) ]
/ dt   \eqno(252)
$$

\vfill
\eject
\vskip 15pt
\vbox{
\hang \noindent {{\it Table 7.} 
History Fields}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Name & Description & Units \cr
\noalign{\hrule height 1pt}
         tmask        & T grid mask (0 = land, 1 = ocean)               &             \cr
         tarea ($A^t$)& area of T grid cells                            & m$^2$       \cr
         uarea ($A^u$)& area of U grid cells                            & m$^2$       \cr
         dxu          & U cell grid width longitudinally through middle & m           \cr
         dyu          & U cell grid width latitudinally through middle  & m           \cr
         dxt          & T cell grid width longitudinally through middle & m           \cr
         dyt          & T cell grid width latitudinally through middle  & m           \cr
         HTN          & T cell width on north side                      & m           \cr
         HTE          & T cell width on east  side                      & m           \cr
         ANGLE  ($\chi^u$)  & angle grid makes with lat line on U grid  & radians     \cr
         ANGLET ($\chi^t$)  & angle grid makes with lat line on T grid  & radians     \cr
         hi ($V$)            & grid box mean ice thickness                   & m              \cr
         hs ($V_s$)          & grid box mean snow thickness                  & m              \cr
         Tsfc ($T_s$)        & snow/ice surface temperature (cpl)            & $^\circ$C      \cr
         aice ($A$)          & aggregate ice area (cpl)                      & \%             \cr
         u ($u_i$)           & x direction ice velocity                      & cm s$^{-1}$    \cr
         v ($v_i$)           & y direction ice velocity                      & cm s$^{-1}$    \cr
         fswdn ($F_{SWDN}$)  & down solar flux                               & W m$^{-2}$     \cr
         flwdn ($F_{LWDN}$)  & down longwave flux                            & W m$^{-2}$     \cr 
         snow ($F_{SNW}$)    & snowfall rate (cpl)                           & cm day$^{-1}$  \cr
    snow$_{ai}$ ($F_{SNW}A$) & snowfall rate                                 & cm day$^{-1}$  \cr
         rain ($F_{RN}$)     & rainfall rate (cpl)                           & cm day$^{-1}$  \cr
    rain$_{ai}$ ($F_{RN}A$)  & rainfall rate                                 & cm day$^{-1}$  \cr
         sst ($T_o$)         & sea surface temperature                       & $^\circ$C      \cr
         sss ($S_o$)         & sea surface salinity                          & ppt            \cr
         uocn ($u_o$)        & x direction ocean current                     & cm s$^{-1}$    \cr
         vocn ($v_o$)        & y direction ocean current                     & cm s$^{-1}$    \cr
         frzmlt ($F_{Qoi}$)  & freezing/melting potential                    & W m$^{-2}$     \cr 
         fswabs ($F_{SW}$)   & snow/ice/ocn absorbed solar flux (cpl)        & W m$^{-2}$     \cr
 fswabs$_{ai}$ ($F_{SW}A$)   & snow/ice/ocn absorbed solar flux              & W m$^{-2}$     \cr
         albsni ($\alpha_{bb}$) & snow-ice broad band albedo                 & \%             \cr
         alvdr  ($\alpha_{vsdr}$) & visible direct albedo                    & \%             \cr
         alidr  ($\alpha_{nidr}$) & near IR direct albedo                    & \%             \cr
\noalign{\hrule height 1pt}
}}}
\vskip 15pt
\vbox{
\hang \noindent {{\it Table 7 continued.} 
History Fields}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Name & Description & Units \cr
\noalign{\hrule height 1pt}
         flat ($F_{LH}$)     & latent heat flux (cpl)                        & W m$^{-2}$     \cr
  flat$_{ai}$ ($F_{LH}A$)    & latent heat flux                              & W m$^{-2}$     \cr
         fsens ($F_{SH}$)    & sensible heat flux (cpl)                      & W m$^{-2}$     \cr
 fsens$_{ai}$  ($F_{SH}A$)   & sensible heat flux                            & W m$^{-2}$     \cr
         flwup ($F_{LWUP}$)  & upward longwave flux (cpl)                    & W m$^{-2}$     \cr
 flwup$_{ai}$ ($F_{LWUP}A$)  & upward longwave flux                          & W m$^{-2}$     \cr
         evap ($F_{EVAP}$)   & evaporative water flux (cpl)                  & cm day$^{-1}$  \cr
  evap$_{ai}$ ($F_{EVAP}A$)  & evaporative water flux                        & cm day$^{-1}$  \cr
         Tref ($T_{REF}$)    & 2m reference temperature (cpl)                & $^\circ$C      \cr
         Qref ($T_{REF}$)    & 2m reference specific humidity (cpl)          & g/kg           \cr
congel ($\Sigma (\delta h\vert_{\rm basal}>0)_n A_n$) & congelation ice growth & cm day$^{-1}$\cr
frazil ($F_{Qoi}/q_f >0$)    & frazil ice growth                             & cm day$^{-1}$  \cr
snoice ($\Sigma \vert z_{int}\vert_n A_n$) & snow-ice conversion             & cm day$^{-1}$  \cr 
meltt  ($\Sigma (\delta h\vert_{\rm melt}<0)_n A_n$)  & top ice melt         & cm day$^{-1}$  \cr
meltb  ($\Sigma (\delta h\vert_{\rm basal}<0)_n A_n$) & basal ice melt       & cm day$^{-1}$  \cr
meltl  ($VR_{side}$)         & lateral ice melt                              & cm day$^{-1}$  \cr
         fresh ($F_{Wo}$)    & freshwater flux ice to ocean (cpl)            & cm day$^{-1}$  \cr
   fresh$_{ai}$ ($F_{Wo}A$)  & freshwater flux ice to ocean                  & cm day$^{-1}$  \cr
         fsalt ($F_{So}$)    & salt flux ice to ocean (cpl)         & kg m$^{-2}$ day$^{-1}$  \cr
   fsalt$_{ai}$ ($F_{So}A$)  & salt flux ice to ocean               & kg m$^{-2}$ day$^{-1}$  \cr
         fhnet ($F_{Qio}$)   & heat flux ice to ocean (cpl)                  & W m$^{-2}$     \cr
   fhnet$_{ai}$ ($F_{Qio}A$) & heat flux ice to ocean                        & W m$^{-2}$     \cr
       fswthru ($F_{SWo}$)   & SW thru ice to ocean (cpl)                    & W m$^{-2}$     \cr
 fswthru$_{ai}$ ($F_{SWo}A$) & SW thru ice to ocean                          & W m$^{-2}$     \cr
         strairx ($\tau_{ax}$)  & x direction atm/ice stress (cpl)           & N m$^{-2}$     \cr 
         strairy ($\tau_{ay}$)  & y direction atm/ice stress (cpl)           & N m$^{-2}$     \cr
  strtltx ($H_{ox}={\bar m}g\partial H_o / \partial x$)& x direction sea surface tilt stress & N m$^{-2}$ \cr
  strtlty ($H_{oy}={\bar m}g\partial H_o / \partial y$)& y direction sea surface tilt stress & N m$^{-2}$ \cr
  strcorx ($+{\bar m}fv_i$)      & x direction coriolis stress               & N m$^{-2}$     \cr
  strcory ($-{\bar m}fu_i$)      & y direction coriolis stress               & N m$^{-2}$     \cr
         strocnx ($\tau_{ox}$)   & x direction ocean/ice stress (cpl)        & N m$^{-2}$     \cr
         strocny ($\tau_{oy}$)   & y direction ocean/ice stress (cpl)        & N m$^{-2}$     \cr
\noalign{\hrule height 1pt}
}}}
\vskip 15pt
\vskip 15pt
\vbox{
\hang \noindent {{\it Table 7 continued.} 
History Fields}
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Name & Description & Units \cr
\noalign{\hrule height 1pt}
strintx $(\nabla\cdot\sigma)_x$ & x direction div internal ice stress tensor & N m$^{-2}$ \cr 
strinty $(\nabla\cdot\sigma)_y$ & y direction div internal ice stress tensor & N m$^{-2}$ \cr 
         strength ($P$)        & compressive ice strength                    & N m$^{-1}$     \cr
         divu ($\dot\epsilon_I$)& strain rate (divergence)                   & \% day$^{-1}$  \cr
         shear ($\dot\epsilon_{II}$)& strain rate (shear)                    & \% day$^{-1}$  \cr
       sig1 ($\sigma_I /P$)  & normalized principal stress component 1       &                \cr
       sig2 ($\sigma_{II}/P$)& normalized principal stress component 2       &                \cr 
         dvidtt $(\partial V / \partial t)_T$ & ice volume tendency thermodynamics & cm day$^{-1}$\cr
         dvidtd $(\partial V / \partial t)_D$ & ice volume tendency dynamics   & cm day$^{-1}$  \cr
         daidtt $(\partial A / \partial t)_T$ & area tendency thermodynamics   & \% day$^{-1}$  \cr
         daidtd $(\partial A / \partial t)_D$ & area tendency dynamics         & \% day$^{-1}$  \cr
opening $({\dot\epsilon_I}+\vert\dot\epsilon\vert C(\theta))$ & lead opening rate & \% day$^{-1}$ \cr
         aice1 ($A_1$)       & ice area (category 1)                         & \%             \cr
         aice2 ($A_2$)       & ice area (category 2)                         & \%             \cr
         aice3 ($A_3$)       & ice area (category 3)                         & \%             \cr
         aice4 ($A_4$)       & ice area (cateogry 4)                         & \%             \cr
         aice5 ($A_5$)       & ice area (category 5)                         & \%             \cr
\noalign{\hrule height 1pt}
}}}
\vskip 15pt

\vfill
\eject
\noindent {\bf 7. Summary}

During 1999-2003, the CCSM Polar Climate Working Group made several
recommendation for improvement to the sea ice model. Additionally,
during 2003-2004 the broader CCSM community imposed requirements 
concerning numerical efficiency and code architecture. The resulting 
CCSM3 sea ice model, CSIM5, addresses these recommendations and
requirements, and includes:

\parskip 10pt
\parindent 15pt

(1) a plastic rheology with an elliptical yield curve

(2) enhanced sea ice thermodynamics 

(3) an ice thickness distribution

(4) elimination of spurious polar convergence near the north pole

(5) an ice model on same grid as the ocean model

(6) an efficient parallel and vector architecture

(7) an active ice only framework for testing 

\parindent 0pt

CSIM5 represents a significant improvement to the original 
sea ice model of CSM1.

\vskip 15pt
\noindent
{\bf Acknowledgments}

Thanks are due to the numerous individals, both at NCAR and in the
scientific community, who contributed to the development of the
CCSM3 sea ice component, CSIM5.

\vfill
\eject
\vskip 15pt
\noindent

{\bf List of Physical Constants}
\vskip 15pt
\vbox{
\hang \noindent { }
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule & 
            \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule & 
            \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Symbol & Code Symbol   &  Description & Value \cr
\noalign{\hrule height 1pt}
         $\rho_s$  &rhos       & density of snow                    & 330 kg m$^{-3}$           \cr
         $\rho_i$  &rhoi       & density of ice                     & 917 kg m$^{-3}$           \cr
         $\rho_o$  &rhow       & density of seawater                & 1026 kg m$^{-3}$          \cr 
         $C_p$     &cp\_air    & specific heat of atmosphere dry    & 1005 J kg$^{-1}$ K$^{-1}$ \cr
         $C_{pwv}$ &cp\_wv     & specific heat of atmosphere water  & 1810 J kg$^{-1}$ K$^{-1}$ \cr
         $c_s$     &cp\_sno    & specific heat of snow              & 0    J kg$^{-1}$ K$^{-1}$ \cr
         $c_0$     &cp\_ice    & specific heat of fresh ice         & 2054 J kg$^{-1}$ K$^{-1}$ \cr
         $c_o$     &cp\_ocn    & specific heat of ocean             & 4218 J kg$^{-1}$ K$^{-1}$ \cr
         $z_i$     &iceruf     & aerodynamic roughness of ice       & 5.0x10$^{-4}$ m \cr
         $z_{ref}$ &zref       & reference height for bulk fluxes   & 10 m     \cr
         $q_1(ice)$&qqqice     & saturation specific humidity const & 11637800 \cr
         $q_2(ice)$&TTTice     & saturation specific humidity const & 5897.8   \cr
         $q_1(ocean)$&qqqocn   & saturation specific humidity const & 627572.4 \cr
         $q_2(ocean)$&TTTocn   & saturation specific humidity const & 5107.4   \cr
         $c_d$     &dragw      & drag coefficient for water on ice  & 0.00536  \cr
         $h_{min}$ &hi\_min    & minimum ice thickness for cat 1    & 0.1 m \cr
         $h_{smin}$&hs\_min    & minimum snow depth for heat eqn    & 0.01 m or 0.00001 m \cr
         $k_{s}$   &ksno       & thermal conductivity of snow       & 0.31  W m$^{-1}$ K$^{-1}$\cr
         $k_{fi}$  &kice       & thermal conductivity of fresh ice  & 2.0340 W m$^{-1}$ K$^{-1}$\cr
         $\beta$   &betak      & thermal conductivity ice constant  & 0.1172 W m$^{-1}$ ppt$^{-1}$ \cr
         $u_{min}$ &ustar\_min & minimum ice/ocean friction velocity& 0.001ms$^{-1}$ \cr
         $L_{i}$   &Lfresh    & latent heat of fusion of ice       & 3.337x10$^5$ J kg$^{-1}$ \cr
         $L_{s}$   &Lsub      & latent heat of sublimation         & 2.835x10$^6$ J kg$^{-1}$ \cr
         $L_{v}$   &Lvap      & latent heat of vaporization        & 2.501x10$^6$ J kg$^{-1}$ \cr
         $T_{f}$   &Tffresh   & freezing temperature of freshwater & 273.15 K          \cr
         $T_{of}$  &Tf        & freezing temperature of ocean      & -1.8$^\circ$C       \cr
         $T_{melt}$&Timelt,Tsmelt & melting temperature of top surface & 0 $^\circ$C          \cr
         $\Delta T_{err}$&Tsf\_errmax & maximum error tolerance            & $5 \times 10^{-4}$ $^\circ$C  \cr
         $\mu$     &depressT  & ocean freezing temperature constant& 0.054 $^\circ$C ppt$^{-1}$ \cr
\noalign{\hrule height 1pt}
}}}
\vskip 15pt
\eject

{\bf Constants continued}
\vskip 15pt
\vbox{
\vskip .5true cm
{\offinterlineskip
\halign{\vrule width 1.5pt height 10pt depth 8pt
            \ \noindent # \hfil \vrule
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule & 
            \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
          & \ \hfil # \hfil \vrule & \ \hfil # \hfil \vrule & 
            \ \hfil # \hfil \vrule width 1.5pt height 10pt depth 8pt
            width 1.5pt \cr
\noalign{\hrule height 1pt}
 Symbol & Code Symbol   &  Description & Value \cr
\noalign{\hrule height 1pt}
         $S_o$     &ocn\_ref\_salinity  & ocean frazil ice ref salinity      & 34.7 ppt    \cr
         $S_i$     &ice\_ref\_salinity  & sea ice frazil ice ref salinity    & 4 ppt       \cr
         $g_e$     &gravit    & gravitational acceleration         & 9.80616  m s$^{-2}$           \cr      
         $\sigma_{sb}$&stefan\_boltzmann   & Stefan-Boltzmann constant & 5.67x10$^{-8}$ W m$^{-2}$ K$^{-4}$   \cr 
         $\varepsilon$&emissivity  & snow/ice emissivity           & 0.95           \cr
         $\kappa_{vs}$&kappav & ice SW visible extinction coefficient & 1.4 m$^{-1}$   \cr
         $\kappa_{ni}$&kappan & ice SW near-ir extinction coefficient & 17.6 m$^{-1}$  \cr 
         $\alpha_o$   &albocn & ocean albedo                       & 0.06  \cr 
         $E_0$        &eyc    & dynamic constant                   & 0.36  \cr 
         $C_s$        &---    & fraction of shear used for ridging & 0.225 \cr 
         $G^\ast$     &gstar  & accumulative ice fraction for ridging & 0.15 \cr 
         $K$          &cK     & max ridged ice thickness constant  & 100 m \cr 
         $Z$          &Zfric  & ratio E dissipated to potential E gain     & 17    \cr 
         $C_{pe}$     &cpe    & ice/ocn potential energy constant  & 450 N m$^{-3}$ \cr
\noalign{\hrule height 1pt}
}}}
\vskip 15pt

\vfill
\eject
\noindent
{\bf Appendix}

During the development of CSIM5, several physics options were incorporated into the 
sea ice model. Some of these relate to options not chosen as default, and others 
are possible candidates for future versions of CSIM. These physics options have been 
retained for the released version of CSIM5. In this appendix these options are briefly
described.

It is useful to run CSIM5 in column mode, namely, as
a {\bf 1D column model} with prescribed inputs (i.e. SHEBA data or other). This 
allows efficient testing and validation of the vertical sea ice physics.
This option has not been thoroughly tested and is not supported.

An earlier method for evaluating thickness space transport is the {\bf delta scheme}
of Bitz (2000) and Bitz et al. (2001). This method represents the distribution function
$g(h)$ as a sum of delta functions, one for each populated ice thickness 
category. This physics option has the limitation of underpopulating the ITD, 
resulting in jumps in properties across category boundaries. The linear remapping
scheme in CSIM5 was chosen as default because it represents $g(h)$ in each category
as a linear function, and allows for incremental transport between categories, yielding
a much smoother representation.

The previous version of the sea ice model evaluated the ice strength using the
method of Hibler (1979). This {\bf ice strength parameterization} depends only
on the mean ice thickness and ice fraction. The scheme in CSIM5 makes use of 
the ice participating in ridging, but is dependable only if the ITD is well 
resolved (i.e. has at least five categories). For fewer categories,
it is preferable to use the Hibler (1979) expression.

Very high spatial resolution implementations of the EVP dynamics using a reasonable 
subcycling time step can produce simulations for which the elastic waves useful for 
regularization are not sufficiently damped. In this case, a physical option of
{\bf damping of elastic waves} is available, as in Hunke (2001) and
Hunke and Dukowicz (2002). The
elastic wave damping is enhanced by reducing the ice strength compared to the
default in CSIM5. For the resolution and subcycling time steps used in CSIM5, 
this option is not necessary.

The CSIM5 ridging scheme by default places snow on ice participating in ridging
into the ocean. If {\bf snow into ocean} is set false, then snow cover is 
retained on ridged ice.

The previous horizontal advection scheme was MPDATA, and before that upwind.
Both {\bf upwind advection} and {\bf MPDATA} are retained as options in CSIM5. 
These are useful for assessing differences among advection schemes.

\vfill\eject

\vskip 10pt
\centerline {\bf References}

\def\refpar{\par\hangindent=3em\hangafter=1}
\def\reference{\relax\refpar}

\reference
Allison, I., R.~E. Brandt, and S.~G. Warren, 1993: East Antarctic sea ice:
albedo, thickness distribution, and snow cover. 
  {\it J. Geophys. Research.}, {\bf 98}, 12417--12429.

\reference
Bettge, T.~W., J.~W.Weatherly, W.~M.Washington, D.Pollard, B.~P.Briegleb,
W.~G.Strand Jr., 1996: The NCAR CSM Sea Ice Model. NCAR/TN-425+STR pp 25.

\reference
Bitz, C.~M., 2000: Documentation of a Lagrangian sea ice thickness distribution
  model with energy-conserving thermodyanmics. U. of Washington  APL-UW TM 4-99 

\reference
Bitz, C.~M., M.~Holland, M.~Eby and A.~J. Weaver, 2001: Simulating the
  ice-thickness distribution in a coupled climate model. 
  {\it J. Geophys. Res.}, {\bf 106}, 2441--2464.

\reference
Bitz, C.~M. and W.~H. Lipscomb, 1999: An energy-conserving thermodynamic model
  of sea ice. {\it J. Geophys. Res.}, {\bf 104}, 15,669--15,677.

\reference
Boville, B.~A. and P.~R. Gent, 1998: The NCAR Climate System Model, Version One 
  {\it J. Climate}, {\bf 11}, 1115-1130.

\reference
Briegleb, B.~P., C.~M. Bitz, E.~C. Hunke, W.~H. Lipscomb, and J.~L. Schramm, 2002:
  Description of the Community Climate System Model Version 2 Sea Ice Model,
  CCSM Web Page- http://www.ccsm.ucar.edu/models/ccsm2.0.1/csim/

\reference
Connolley, W.~M., J.~M. Gregory, E.~C. Hunke, and A.~J. McLaren, 2004: On the
  consistent scaling of terms in the sea ice dynamics equation. {\it J. Phys. Oceanogr.}, 
  in press. LA-UR-03-4772.

\reference
Curry, J.~A., J.~L. Schramm, E.~E. Ebert, 1995: Sea ice-albedo feedback mechanism.
  {\it J. Climate}, {\bf 8}, 240--247.

\reference
Curry, J.~A., J.~L. Schramm, D.~K. Perovich, and J.~O. Pinto, 2001: Applications of
SHEBA/FIRE data to evaluation of snow/ice albedo parameterizations.
  {\it J. Geophys. Research.}, {\bf 106}, D14, 15345-15355.

\reference
Dukowicz, J.~K. and J.~R. Baumgardner, 2000: Incremental remapping as a 
transport/advection algorithm.
  {\it J. Comput. Phys.}, {\bf 160}, 318-335.

\reference
Ebert, E.~E. and J.~A. Curry, 1993: An intermediate one-dimensional sea ice
model for investigating ice-atmosphere interactions.
  {\it J. Geophys. Research.}, {\bf 98}, 10085--10109.

\reference
Flato, G.~M. and W.~D. Hibler, III, 1992: Modeling pack ice as a cavitating fluid.
  {\it J. Phys. Oceanogr.}, {\bf 22}, 626--651.

\reference
Flato, G.~M. and W.~D. Hibler, III, 1995: Ridging and strength in modelling the
  thickness distribution of Arctic sea ice.
  {\it J. Geophys. Research.}, {\bf C9}, 18611--18626.

\reference
Grenfell, T.~C., S.~G. Warren, and P.~C. Mullen, 1994: Reflection of solar radiation
by the Antarctic snow surface at ultraviolet, visible, and near-infrared wavelengths.
  {\it J. Geophys. Research.}, {\bf 99}, 18669--18684.

\reference
Hibler, W.~D.,III, 1979: A dynamic thermodynamic sea ice model.
  {\it J. Phys. Oceanogr.}, {\bf 9}, 815--846.

\reference
Hibler, W.~D.,III, 1980: Modeling a variable thickness ice cover.
  {\it Mon. Wea. Rev.}, {\bf 108}, 1943--1973.

\reference
Hogstrom, U., 1988: Non-dimensional wind and temperature profiles in the
atmospheric surface layer: a re-evaluation. 
  {\it Boundary-Layer Meteorology}, {\bf 42}, 55--78.

\reference
Hopkins, M.~A. and W.~D. Hibler, III, 1991: On the ridging of a thin sheet of lead ice.
  {\it Annals of Glaciology}, {\bf 15}, 81--86.

\reference
Hunke, E.~C. and J.~K. Dukowicz, 1997: An elastic-viscous-plastic model for sea
  ice dynamics. {\it J. Phys. Oceanogr.}, {\bf 27}, 1849--1867.

\reference
Hunke, E.~C. and W.~H. Lipscomb, 2004: CICE: the Los Alamos sea ice model, documentation and software
  User's Manual. T-3 Fluid Dynamics Group, Los Alamos National Laboratory, Tech. Rep. LA-CC-98-16 v.3.1

\reference
Hunke, E.~C., 2001: Viscous-plastic sea ice dynamics with the evp model: linearization issues.
  {\it J. Comp. Phys.}, {\bf 170}, 18--38.

\reference
Hunke, E.~C. and J.~K. Dukowicz, 2002: The elastic-viscous-plastic
  sea ice dynamics model in general orthogonal curvilinear coordinates
  on a sphere---incorporation of metric terms. {\it Monthly Weather Review},
  {\bf 130}, 1848--1865.

\reference
Hunke, E.~C. and J.~K. Dukowicz, 2003: The sea ice momentum equation in the free drift regime.
 Tech. Rep. LA-UR-03-2219, Los Alamos National Laboratory.

\reference
Kreyscher, M., M.Harder, P.Lemke and G.M.Flato., 2000: Results of the
  Sea Ice Model Intercomparison Project. {\it J. Geophys. Research.}, {\bf 105} 
  C5, 11299--11320.

\reference
Large, W.G., 1998: Modeling and parameterizing the ocean planetary boundary layer. 
{\it Ocean Modeling and Parameterization}, 81--120, E.P.Chassignet and J.Verron (eds.),
Kluwer Academic Publishers. Printed in the Netherlands.

\reference
Lipscomb, W.~H, 2001: Remapping the thickness distribution in sea ice models.
  {\it J. Geophys. Research.}, {\bf 106}, 13,989--14,000.

\reference
Lipscomb, W.~H and E.~C. Hunke, 2004: Modeling sea ice transport using
incremental remapping. {\it Mon. Wea. Rev.}, {\bf 132}, 1341--1354.

\reference
Maykut, G.~A. and N. Untersteiner, 1971: Some results from a time-dependent
thermodynamic model of sea ice. 
  {\it J. Geophys. Research.}, {\bf 76}, 1550--1575.

\reference
Maykut, G.~A. and D. Perovich, 1987: The role of shortwave radiation in the 
summer decay of sea ice cover.
  {\it J. Geophys. Research.}, {\bf 92}, 7032--7044.

\reference
McPhee, M.~G., 1992: Turbulent heat flux in the upper ocean under sea ice. 
  {\it J. Geophys. Research.}, {\bf 97}, 5365-5379.

\reference
The NCAR CSM Flux Coupler, 1996. NCAR/TN-424+STR 

\reference
Ono, M., 1967: Specific heat and heat of fusion of sea ice, in {\it Physics
  of Snow and Ice}, edited by H. Oura, {\bf 1}, 599-610, 
  Inst. of Low Temp. Sci., Hokkaido, Japan.

\reference
Paulson, C.~A. and J.~J. Simpson, 1977: Irradiance measurements in the upper ocean.
  {\it J. Phys. Oceanogr.}, {\bf 7}, 952--956.

\reference
Rothrock, D.~A., 1975: The energetics of the plastic deformation of pack ice by ridging.
  {\it J. Geophys. Research.}, {\bf 80}, 4514-4519.

\reference
Rothrock, D.~A. and A.~S. Thorndike, 1984: Measuring the sea ice floe size distribution.
  {\it J. Geophys. Research.}, {\bf 89}, 6477--6486.

\reference
Semtner, A.~J., 1976: A model for the thermodynamic growth of sea ice in
  numerical investigations of climate. {\it J. Phys. Oceanogr.}, {\bf 6},
  379--389.

\reference
Smolarkiewicz, P.~K., 1984: A fully multidimensional positive definite
  advection scheme with small implicit diffusion. {\it J. Comput. Phys.},
  {\bf 54 }, 325--362.

\reference
Stern, H.~L.,D.~A. Rothrock, and R. Kwok, 1995: Open water production in 
  Arctic sea ice: satellite measurements and model parameterizations.
  {\it J. Geophys. Res.}, {\bf 100}, 20601--20612.

\reference
Thorndike, A.~S., D.~A. Rothrock, G.~A. Maykut and R.~Colony, 1975: The
  thickness distribution of sea ice. {\it J. Geophys. Res.}, {\bf 80},
  4501--4513.

\reference
Untersteiner, N., 1961: On the mass and heat budget of arctic sea ice.
  {\it Arch. Meteorol. Geophys. Bioklimatol.,A} {\bf 12}, 151--182.

\reference
Weatherly, J.~W., B.~P. Briegleb, W.~G. Large and J.~A. Maslanik, 1998: 
  Sea Ice and Polar Climate in the NCAR CSM.
  {\it J. Climate}, {\bf 11}, 1472--1486.

\vfill\eject

\bye
